DTU on simulated datasets
Kimon Froussios
03/06/2019
Origin of simulated datasets
Soneson et al., Genom Biol 2016 17:12 : https://doi.org/10.1186/s13059-015-0862-3 ArrayExpress repository of the simulated datasets: https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3766/
Prelims
# libraries
library(rats)
library(DRIMSeq)
library(data.table)
library(ggplot2)
library(matrixStats)
library(parallel)
library(gganimate)
library(dplyr)
# ggplot2 themes
gth <- theme(axis.line.x= element_line(),
axis.line.y= element_line(),
strip.background= element_rect(fill= "grey95"),
strip.text.y= element_text(size= rel(1.2)),
strip.text.x= element_text(size= rel(1.1)),
panel.grid.major.x= element_blank(),
panel.grid.minor.x= element_blank(),
panel.grid.major.y= element_line(colour='grey90'),
panel.grid.minor.y= element_blank(),
panel.background= element_rect(fill = "white"),
legend.key = element_rect(fill = 'white'),
panel.spacing = unit(1, "lines") )
gth2 <- theme(axis.line.x= element_line(),
axis.text.x=element_text(angle=90),
axis.line.y= element_line(),
strip.background= element_rect(fill= "grey95"),
strip.text.y= element_text(size= rel(1.2)),
strip.text.x= element_text(size= rel(1.1)),
panel.grid.major.x= element_line(colour='grey90'),
panel.grid.minor.x= element_line(colour='grey95'),
panel.grid.major.y= element_line(colour='grey90'),
panel.grid.minor.y= element_line(colour='grey95'),
panel.background= element_rect(fill = "white"),
legend.key = element_rect(fill = 'white'),
panel.spacing = unit(1, "lines") )
gth3 <- theme(axis.line.x= element_line(),
axis.text.x=element_text(angle=90, hjust=1),
axis.line.y= element_line(),
strip.background= element_rect(fill= "grey95"),
strip.text.y= element_text(size= rel(1.2)),
strip.text.x= element_text(size= rel(1.1)),
panel.grid.major.x= element_line(colour='grey95'),
panel.grid.minor.x= element_blank(),
panel.grid.major.y= element_line(colour='grey95'),
panel.grid.minor.y= element_line(colour='grey95'),
panel.background= element_rect(fill = "white"),
legend.key = element_rect(fill = 'white'),
panel.spacing = unit(1, "lines"),
strip.placement = "outside")Import data
basedir <- '~'
# Import SUPPA2 results
sdk <- fread(file.path(basedir, 'DE', 'suppa_soneson_Dm70-kallisto.tsv.dpsi.temp.0'))
sds <- fread(file.path(basedir, 'DE', 'suppa_soneson_Dm70-salmon.tsv.dpsi.temp.0'))
shk <- fread(file.path(basedir, 'DE', 'suppa_soneson_Hs71-kallisto.tsv.dpsi.temp.0'))
shs <- fread(file.path(basedir, 'DE', 'suppa_soneson_Hs71-salmon.tsv.dpsi.temp.0'))
# Import DRIMSeq results
# Gene-level (default)
ddk <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Dm70-kallisto.results.tsv'))
dds <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Dm70-salmon.results.tsv'))
dhk <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Hs71-kallisto.results.tsv'))
dhs <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Hs71-salmon.results.tsv'))
# Transcript-level
ddkt <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Dm70-kallisto.results-t.tsv'))
ddst <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Dm70-salmon.results-t.tsv'))
dhkt <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Hs71-kallisto.results-t.tsv'))
dhst <- fread(file.path(basedir, 'DE', 'drimseq0_soneson_Hs71-salmon.results-t.tsv'))
# Pre-comparison object (contains the abundances)
ddkr <- readRDS(file.path(basedir, 'DE', 'drimseq0_soneson_Dm70-kallisto.RDS'))
ddsr <- readRDS(file.path(basedir, 'DE', 'drimseq0_soneson_Dm70-salmon.RDS'))
dhkr <- readRDS(file.path(basedir, 'DE', 'drimseq0_soneson_Hs71-kallisto.RDS'))
dhsr <- readRDS(file.path(basedir, 'DE', 'drimseq0_soneson_Hs71-salmon.RDS'))
# Import DEXSeq results.
xds <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq0_soneson_Dm70-salmon.RDS')))## Loading required package: DEXSeq## Loading required package: BiocParallel## Loading required package: Biobase## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind,
## colMeans, colnames, colSums, dirname, do.call, duplicated,
## eval, evalq, Filter, Find, get, grep, grepl, intersect,
## is.unsorted, lapply, lengths, Map, mapply, match, mget, order,
## paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind,
## Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which, which.max,
## which.min## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians## The following object is masked from 'package:DRIMSeq':
##
## samples## Loading required package: SummarizedExperiment## Loading required package: GenomicRanges## Loading required package: stats4## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:dplyr':
##
## first, rename## The following objects are masked from 'package:data.table':
##
## first, second## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice## The following object is masked from 'package:data.table':
##
## shift## Loading required package: GenomeInfoDb## Loading required package: DelayedArray##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following objects are masked from 'package:base':
##
## aperm, apply## Loading required package: DESeq2##
## Attaching package: 'DESeq2'## The following object is masked from 'package:DRIMSeq':
##
## results## Loading required package: AnnotationDbi##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: RColorBrewerxhs <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq0_soneson_Hs71-salmon.RDS')))
xsds <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq0-stageR_soneson_Dm70-salmon.RDS')))
xshs <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq0-stageR_soneson_Hs71-salmon.RDS')))
# Size-normalised abundances
xdsc <- as.data.table(counts(readRDS(file.path(basedir, 'DE', 'dexseq0_soneson_Dm70-salmon.RDS')), normalized=TRUE))
xhsc <- as.data.table(counts(readRDS(file.path(basedir, 'DE', 'dexseq0_soneson_Hs71-salmon.RDS')), normalized=TRUE))
# Import ground truth
dtruth <- fread(file.path(basedir, 'truth_drosophila_non_null_missing20_ms.txt'))
htruth <- fread(file.path(basedir, 'truth_human_non_null_missing20_ms.txt'))
# Import RATs results with 0.6.4 default thresholds
rdk <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Dm70-kallisto_def.RDS'))
rds <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Dm70-salmon_def.RDS'))
rhk <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Hs71-kallisto_def.RDS'))
rhs <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Hs71-salmon_def.RDS'))
# Import RATs results with different thresholds.
rdkth <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Dm70-kallisto_th4.RDS'))
rdsth <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Dm70-salmon_th4.RDS'))
rhkth <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Hs71-kallisto_th4.RDS'))
rhsth <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Hs71-salmon_th4.RDS'))
rdkth2 <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Dm70-kallisto_th5.RDS'))
rdsth2 <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Dm70-salmon_th5.RDS'))
rhkth2 <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Hs71-kallisto_th5.RDS'))
rhsth2 <- readRDS(file.path(basedir, 'DE', 'rats_soneson_Hs71-salmon_th5.RDS'))# Import SUPPA2 results
sdk10 <- fread(file.path(basedir, 'DE', 'suppa10_soneson_Dm70-kallisto.tsv.dpsi.temp.0'))
sds10 <- fread(file.path(basedir, 'DE', 'suppa10_soneson_Dm70-salmon.tsv.dpsi.temp.0'))
shk10 <- fread(file.path(basedir, 'DE', 'suppa10_soneson_Hs71-kallisto.tsv.dpsi.temp.0'))
shs10 <- fread(file.path(basedir, 'DE', 'suppa10_soneson_Hs71-salmon.tsv.dpsi.temp.0'))
# Import DRIMSeq results
# Gene-level (default)
ddk10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Dm70-kallisto.results.tsv'))
dds10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Dm70-salmon.results.tsv'))
dhk10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Hs71-kallisto.results.tsv'))
dhs10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Hs71-salmon.results.tsv'))
# Transcript-level
ddkt10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Dm70-kallisto.results-t.tsv'))
ddst10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Dm70-salmon.results-t.tsv'))
dhkt10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Hs71-kallisto.results-t.tsv'))
dhst10 <- fread(file.path(basedir, 'DE', 'drimseq10_soneson_Hs71-salmon.results-t.tsv'))
# Pre-comparison object (contains the abundances)
ddkr10 <- readRDS(file.path(basedir, 'DE', 'drimseq10_soneson_Dm70-kallisto.RDS'))
ddsr10 <- readRDS(file.path(basedir, 'DE', 'drimseq10_soneson_Dm70-salmon.RDS'))
dhkr10 <- readRDS(file.path(basedir, 'DE', 'drimseq10_soneson_Hs71-kallisto.RDS'))
dhsr10 <- readRDS(file.path(basedir, 'DE', 'drimseq10_soneson_Hs71-salmon.RDS'))
# # Import DEXSeq results.
xds10 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq10_soneson_Dm70-salmon.RDS')))
xhs10 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq10_soneson_Hs71-salmon.RDS')))
xsds10 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq10-stageR_soneson_Dm70-salmon.RDS')))
xshs10 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq10-stageR_soneson_Hs71-salmon.RDS')))
# Size-normalised abundances
xdsc10 <- as.data.table(counts(readRDS(file.path(basedir, 'DE', 'dexseq10_soneson_Dm70-salmon.RDS')), normalized=TRUE))
xhsc10 <- as.data.table(counts(readRDS(file.path(basedir, 'DE', 'dexseq10_soneson_Hs71-salmon.RDS')), normalized=TRUE))
# Import RATs results with 0.6.4 default thresholds
rdk10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Dm70-kallisto_def.RDS'))
rds10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Dm70-salmon_def.RDS'))
rhk10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Hs71-kallisto_def.RDS'))
rhs10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Hs71-salmon_def.RDS'))
# Import RATs results with different thresholds.
rdkth10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Dm70-kallisto_th4.RDS'))
rdsth10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Dm70-salmon_th4.RDS'))
rhkth10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Hs71-kallisto_th4.RDS'))
rhsth10 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Hs71-salmon_th4.RDS'))
rdkth210 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Dm70-kallisto_th5.RDS'))
rdsth210 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Dm70-salmon_th5.RDS'))
rhkth210 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Hs71-kallisto_th5.RDS'))
rhsth210 <- readRDS(file.path(basedir, 'DE', 'rats10_soneson_Hs71-salmon_th5.RDS'))# Import SUPPA2 results
sdk50 <- fread(file.path(basedir, 'DE', 'suppa50_soneson_Dm70-kallisto.tsv.dpsi.temp.0'))
sds50 <- fread(file.path(basedir, 'DE', 'suppa50_soneson_Dm70-salmon.tsv.dpsi.temp.0'))
shk50 <- fread(file.path(basedir, 'DE', 'suppa50_soneson_Hs71-kallisto.tsv.dpsi.temp.0'))
shs50 <- fread(file.path(basedir, 'DE', 'suppa50_soneson_Hs71-salmon.tsv.dpsi.temp.0'))
# Import DRIMSeq results
# Gene-level (default)
ddk50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Dm70-kallisto.results.tsv'))
dds50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Dm70-salmon.results.tsv'))
dhk50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Hs71-kallisto.results.tsv'))
dhs50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Hs71-salmon.results.tsv'))
# Transcript-level
ddkt50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Dm70-kallisto.results-t.tsv'))
ddst50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Dm70-salmon.results-t.tsv'))
dhkt50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Hs71-kallisto.results-t.tsv'))
dhst50 <- fread(file.path(basedir, 'DE', 'drimseq50_soneson_Hs71-salmon.results-t.tsv'))
# Pre-comparison object (contains the abundances)
ddkr50 <- readRDS(file.path(basedir, 'DE', 'drimseq50_soneson_Dm70-kallisto.RDS'))
ddsr50 <- readRDS(file.path(basedir, 'DE', 'drimseq50_soneson_Dm70-salmon.RDS'))
dhkr50 <- readRDS(file.path(basedir, 'DE', 'drimseq50_soneson_Hs71-kallisto.RDS'))
dhsr50 <- readRDS(file.path(basedir, 'DE', 'drimseq50_soneson_Hs71-salmon.RDS'))
# # Import DEXSeq results.
xds50 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq50_soneson_Dm70-salmon.RDS')))
xhs50 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq50_soneson_Hs71-salmon.RDS')))
xsds50 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq50-stageR_soneson_Dm70-salmon.RDS')))
xshs50 <- as.data.table(readRDS(file.path(basedir, 'DE', 'dexseq50-stageR_soneson_Hs71-salmon.RDS')))
# Size-normalised abundances
xdsc50 <- as.data.table(counts(readRDS(file.path(basedir, 'DE', 'dexseq50_soneson_Dm70-salmon.RDS')), normalized=TRUE))
xhsc50 <- as.data.table(counts(readRDS(file.path(basedir, 'DE', 'dexseq50_soneson_Hs71-salmon.RDS')), normalized=TRUE))
# Import RATs results with 0.6.4 default thresholds
rdk50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Dm70-kallisto_def.RDS'))
rds50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Dm70-salmon_def.RDS'))
rhk50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Hs71-kallisto_def.RDS'))
rhs50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Hs71-salmon_def.RDS'))
# Import RATs results with different thresholds.
rdkth50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Dm70-kallisto_th4.RDS'))
rdsth50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Dm70-salmon_th4.RDS'))
rhkth50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Hs71-kallisto_th4.RDS'))
rhsth50 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Hs71-salmon_th4.RDS'))
rdkth250 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Dm70-kallisto_th5.RDS'))
rdsth250 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Dm70-salmon_th5.RDS'))
rhkth250 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Hs71-kallisto_th5.RDS'))
rhsth250 <- readRDS(file.path(basedir, 'DE', 'rats50_soneson_Hs71-salmon_th5.RDS'))In the plots that follow, the different pre-filter values are shown as different frames of the animation.
Preprocess and structure
# Index DRIMSeq by gene
setkey(ddk, gene_id)
setkey(dds, gene_id)
setkey(dhk, gene_id)
setkey(dhs, gene_id)
setkey(ddkt, gene_id)
setkey(ddst, gene_id)
setkey(dhkt, gene_id)
setkey(dhst, gene_id)
# Calculate difference of proporions for DRIMSeq.
txProp <- function (d, asel=c("s1", "s2", "s3"), bsel=c("s4", "s5", "s6")) {
cnts <- as.data.table(counts(d))
setkey(cnts, gene_id, feature_id)
ctsA <- data.table(gene_id=cnts$gene_id, cnt=rowSums(subset(cnts, select=asel)))
ctsB <- data.table(gene_id=cnts$gene_id, cnt=rowSums(subset(cnts, select=bsel)))
ctsA[, total := sum(cnt), by=gene_id]
ctsB[, total := sum(cnt), by=gene_id]
ctsA[, prop := cnt / total]
ctsB[, prop := cnt / total ]
return(ctsB$prop - ctsA$prop)
}
ddkt$Dprop <- txProp(ddkr)
ddst$Dprop <- txProp(ddsr)
dhkt$Dprop <- txProp(dhkr)
dhst$Dprop <- txProp(dhsr)
rm(ddkr)
rm(ddsr)
rm(dhkr)
rm(dhsr)
ddk$maxDprop <- ddkt[, max(abs(Dprop)), by=gene_id]$V1 ; names(ddk)[length(ddk)] <- "maxDprop"
dds$maxDprop <- ddst[, max(abs(Dprop)), by=gene_id]$V1 ; names(dds)[length(dds)] <- "maxDprop"
dhk$maxDprop <- dhkt[, max(abs(Dprop)), by=gene_id]$V1 ; names(dhk)[length(dhk)] <- "maxDprop"
dhs$maxDprop <- dhst[, max(abs(Dprop)), by=gene_id]$V1 ; names(dhs)[length(dhs)] <- "maxDprop"
# Give shorter names to SUPPA2 columns, for ease of typing
names(sdk) <- c("event","dpsi","pval")
names(sds) <- c("event","dpsi","pval")
names(shk) <- c("event","dpsi","pval")
names(shs) <- c("event","dpsi","pval")
# Split SUPPA2 event IDs, to obtain the gene ID and transcript ID
sdk[, gene_id := sapply(as.character(sdk$event), function (x) strsplit(x, ";")[[1]][1])]
sds[, gene_id := sapply(as.character(sds$event), function (x) strsplit(x, ";")[[1]][1])]
shk[, gene_id := sapply(as.character(shk$event), function (x) strsplit(x, ";")[[1]][1])]
shs[, gene_id := sapply(as.character(shs$event), function (x) strsplit(x, ";")[[1]][1])]
sdk[, transc_id := sapply(as.character(sdk$event), function (x) strsplit(x, ";")[[1]][2])]
sds[, transc_id := sapply(as.character(sds$event), function (x) strsplit(x, ";")[[1]][2])]
shk[, transc_id := sapply(as.character(shk$event), function (x) strsplit(x, ";")[[1]][2])]
shs[, transc_id := sapply(as.character(shs$event), function (x) strsplit(x, ";")[[1]][2])]
# Index SUPPA2 by gene
setkey(sdk, gene_id)
setkey(sds, gene_id)
setkey(shk, gene_id)
setkey(shs, gene_id)
# Rename gene and transcript columns
names(xds)[1:2] <- c('geneID', 'txID')
names(xhs)[1:2] <- c('geneID', 'txID')
# Calculate DEXSeq effect size
xdsc[, geneID := xds$geneID]
xhsc[, geneID := xhs$geneID]
xdsc[, txID := xds$txID]
xhsc[, txID := xhs$txID]
xdsc[, sumA := rowSums(xdsc[, 1:3], na.rm=TRUE) ]
xdsc[, sumB := rowSums(xdsc[, 4:6], na.rm=TRUE) ]
xhsc[, sumA := rowSums(xhsc[, 1:3], na.rm=TRUE) ]
xhsc[, sumB := rowSums(xhsc[, 4:6], na.rm=TRUE) ]
xdsc[, totalA := sum(sumA, na.rm=TRUE), by=geneID]
xdsc[, totalB := sum(sumB, na.rm=TRUE), by=geneID]
xhsc[, totalA := sum(sumA, na.rm=TRUE), by=geneID]
xhsc[, totalB := sum(sumB, na.rm=TRUE), by=geneID]
xdsc[, propA := sumA/totalA]
xdsc[, propB := sumB/totalB]
xhsc[, propA := sumA/totalA]
xhsc[, propB := sumB/totalB]
xdsc[(is.nan(propA)), propA := NA_real_] # Replace NaN with NA.
xdsc[(is.nan(propB)), propB := NA_real_]
xhsc[(is.nan(propA)), propA := NA_real_] # Replace NaN with NA.
xhsc[(is.nan(propB)), propB := NA_real_]
xdsc[, Dprop := propB - propA]
xhsc[, Dprop := propB - propA]
xds[, Dprop := xdsc$Dprop]
xhs[, Dprop := xhsc$Dprop]
# Index DEXSeq by gene
setkey(xds, geneID, txID)
setkey(xhs, geneID, txID)
setkey(xsds, geneID, txID)
setkey(xshs, geneID, txID)
xsds[, Dprop := xdsc$Dprop]
xshs[, Dprop := xhsc$Dprop]# Index DRIMSeq by gene
setkey(ddk10, gene_id)
setkey(dds10, gene_id)
setkey(dhk10, gene_id)
setkey(dhs10, gene_id)
setkey(ddkt10, gene_id)
setkey(ddst10, gene_id)
setkey(dhkt10, gene_id)
setkey(dhst10, gene_id)
ddkt10$Dprop <- txProp(ddkr10)
ddst10$Dprop <- txProp(ddsr10)
dhkt10$Dprop <- txProp(dhkr10)
dhst10$Dprop <- txProp(dhsr10)
rm(ddkr10)
rm(ddsr10)
rm(dhkr10)
rm(dhsr10)
ddk10$maxDprop <- ddkt10[, max(abs(Dprop)), by=gene_id]$V1 ; names(ddk10)[length(ddk10)] <- "maxDprop"
dds10$maxDprop <- ddst10[, max(abs(Dprop)), by=gene_id]$V1 ; names(dds10)[length(dds10)] <- "maxDprop"
dhk10$maxDprop <- dhkt10[, max(abs(Dprop)), by=gene_id]$V1 ; names(dhk10)[length(dhk10)] <- "maxDprop"
dhs10$maxDprop <- dhst10[, max(abs(Dprop)), by=gene_id]$V1 ; names(dhs10)[length(dhs10)] <- "maxDprop"
# Give shorter names to SUPPA2 columns, for ease of typing
names(sdk10) <- c("event","dpsi","pval")
names(sds10) <- c("event","dpsi","pval")
names(shk10) <- c("event","dpsi","pval")
names(shs10) <- c("event","dpsi","pval")
# Split SUPPA2 event IDs, to obtain the gene ID and transcript ID
sdk10[, gene_id := sapply(as.character(sdk10$event), function (x) strsplit(x, ";")[[1]][1])]
sds10[, gene_id := sapply(as.character(sds10$event), function (x) strsplit(x, ";")[[1]][1])]
shk10[, gene_id := sapply(as.character(shk10$event), function (x) strsplit(x, ";")[[1]][1])]
shs10[, gene_id := sapply(as.character(shs10$event), function (x) strsplit(x, ";")[[1]][1])]
sdk10[, transc_id := sapply(as.character(sdk10$event), function (x) strsplit(x, ";")[[1]][2])]
sds10[, transc_id := sapply(as.character(sds10$event), function (x) strsplit(x, ";")[[1]][2])]
shk10[, transc_id := sapply(as.character(shk10$event), function (x) strsplit(x, ";")[[1]][2])]
shs10[, transc_id := sapply(as.character(shs10$event), function (x) strsplit(x, ";")[[1]][2])]
# Index SUPPA2 by gene
setkey(sdk10, gene_id)
setkey(sds10, gene_id)
setkey(shk10, gene_id)
setkey(shs10, gene_id)
# Rename gene and transcript columns
names(xds10)[1:2] <- c('geneID', 'txID')
names(xhs10)[1:2] <- c('geneID', 'txID')
# Calculate DEXSeq effect size
xdsc10[, geneID := xds10$geneID]
xhsc10[, geneID := xhs10$geneID]
xdsc10[, txID := xds10$txID]
xhsc10[, txID := xhs10$txID]
xdsc10[, sumA := rowSums(xdsc10[, 1:3], na.rm=TRUE) ]
xdsc10[, sumB := rowSums(xdsc10[, 4:6], na.rm=TRUE) ]
xhsc10[, sumA := rowSums(xhsc10[, 1:3], na.rm=TRUE) ]
xhsc10[, sumB := rowSums(xhsc10[, 4:6], na.rm=TRUE) ]
xdsc10[, totalA := sum(sumA, na.rm=TRUE), by=geneID]
xdsc10[, totalB := sum(sumB, na.rm=TRUE), by=geneID]
xhsc10[, totalA := sum(sumA, na.rm=TRUE), by=geneID]
xhsc10[, totalB := sum(sumB, na.rm=TRUE), by=geneID]
xdsc10[, propA := sumA/totalA]
xdsc10[, propB := sumB/totalB]
xhsc10[, propA := sumA/totalA]
xhsc10[, propB := sumB/totalB]
xdsc10[(is.nan(propA)), propA := NA_real_] # Replace NaN with NA.
xdsc10[(is.nan(propB)), propB := NA_real_]
xhsc10[(is.nan(propA)), propA := NA_real_] # Replace NaN with NA.
xhsc10[(is.nan(propB)), propB := NA_real_]
xdsc10[, Dprop := propB - propA]
xhsc10[, Dprop := propB - propA]
xds10[, Dprop := xdsc10$Dprop]
xhs10[, Dprop := xhsc10$Dprop]
# Index DEXSeq by gene
setkey(xds10, geneID, txID)
setkey(xhs10, geneID, txID)
setkey(xsds10, geneID, txID)
setkey(xshs10, geneID, txID)
xsds10[, Dprop := xdsc10$Dprop]
xshs10[, Dprop := xhsc10$Dprop]# Index DRIMSeq by gene
setkey(ddk50, gene_id)
setkey(dds50, gene_id)
setkey(dhk50, gene_id)
setkey(dhs50, gene_id)
setkey(ddkt50, gene_id)
setkey(ddst50, gene_id)
setkey(dhkt50, gene_id)
setkey(dhst50, gene_id)
ddkt50$Dprop <- txProp(ddkr50)
ddst50$Dprop <- txProp(ddsr50)
dhkt50$Dprop <- txProp(dhkr50)
dhst50$Dprop <- txProp(dhsr50)
rm(ddkr50)
rm(ddsr50)
rm(dhkr50)
rm(dhsr50)
ddk50$maxDprop <- ddkt50[, max(abs(Dprop)), by=gene_id]$V1 ; names(ddk50)[length(ddk50)] <- "maxDprop"
dds50$maxDprop <- ddst50[, max(abs(Dprop)), by=gene_id]$V1 ; names(dds50)[length(dds50)] <- "maxDprop"
dhk50$maxDprop <- dhkt50[, max(abs(Dprop)), by=gene_id]$V1 ; names(dhk50)[length(dhk50)] <- "maxDprop"
dhs50$maxDprop <- dhst50[, max(abs(Dprop)), by=gene_id]$V1 ; names(dhs50)[length(dhs50)] <- "maxDprop"
# Give shorter names to SUPPA2 columns, for ease of typing
names(sdk50) <- c("event","dpsi","pval")
names(sds50) <- c("event","dpsi","pval")
names(shk50) <- c("event","dpsi","pval")
names(shs50) <- c("event","dpsi","pval")
# Split SUPPA2 event IDs, to obtain the gene ID and transcript ID
sdk50[, gene_id := sapply(as.character(sdk50$event), function (x) strsplit(x, ";")[[1]][1])]
sds50[, gene_id := sapply(as.character(sds50$event), function (x) strsplit(x, ";")[[1]][1])]
shk50[, gene_id := sapply(as.character(shk50$event), function (x) strsplit(x, ";")[[1]][1])]
shs50[, gene_id := sapply(as.character(shs50$event), function (x) strsplit(x, ";")[[1]][1])]
sdk50[, transc_id := sapply(as.character(sdk50$event), function (x) strsplit(x, ";")[[1]][2])]
sds50[, transc_id := sapply(as.character(sds50$event), function (x) strsplit(x, ";")[[1]][2])]
shk50[, transc_id := sapply(as.character(shk50$event), function (x) strsplit(x, ";")[[1]][2])]
shs50[, transc_id := sapply(as.character(shs50$event), function (x) strsplit(x, ";")[[1]][2])]
# Index SUPPA2 by gene
setkey(sdk50, gene_id)
setkey(sds50, gene_id)
setkey(shk50, gene_id)
setkey(shs50, gene_id)
# Rename gene and transcript columns
names(xds50)[1:2] <- c('geneID', 'txID')
names(xhs50)[1:2] <- c('geneID', 'txID')
# Calculate DEXSeq effect size
xdsc50[, geneID := xds50$geneID]
xhsc50[, geneID := xhs50$geneID]
xdsc50[, txID := xds50$txID]
xhsc50[, txID := xhs50$txID]
xdsc50[, sumA := rowSums(xdsc50[, 1:3], na.rm=TRUE) ]
xdsc50[, sumB := rowSums(xdsc50[, 4:6], na.rm=TRUE) ]
xhsc50[, sumA := rowSums(xhsc50[, 1:3], na.rm=TRUE) ]
xhsc50[, sumB := rowSums(xhsc50[, 4:6], na.rm=TRUE) ]
xdsc50[, totalA := sum(sumA, na.rm=TRUE), by=geneID]
xdsc50[, totalB := sum(sumB, na.rm=TRUE), by=geneID]
xhsc50[, totalA := sum(sumA, na.rm=TRUE), by=geneID]
xhsc50[, totalB := sum(sumB, na.rm=TRUE), by=geneID]
xdsc50[, propA := sumA/totalA]
xdsc50[, propB := sumB/totalB]
xhsc50[, propA := sumA/totalA]
xhsc50[, propB := sumB/totalB]
xdsc50[(is.nan(propA)), propA := NA_real_] # Replace NaN with NA.
xdsc50[(is.nan(propB)), propB := NA_real_]
xhsc50[(is.nan(propA)), propA := NA_real_] # Replace NaN with NA.
xhsc50[(is.nan(propB)), propB := NA_real_]
xdsc50[, Dprop := propB - propA]
xhsc50[, Dprop := propB - propA]
xds50[, Dprop := xdsc50$Dprop]
xhs50[, Dprop := xhsc50$Dprop]
# Index DEXSeq by gene
setkey(xds50, geneID, txID)
setkey(xhs50, geneID, txID)
setkey(xsds50, geneID, txID)
setkey(xshs50, geneID, txID)
xsds50[, Dprop := xdsc50$Dprop]
xshs50[, Dprop := xhsc50$Dprop]Overview of truth
# Verify which field defines the 1000 transcripts that were intentionally altered
cat("Values of Status")## Values of Statusprint( dtruth[, unique(ds_status)] )## [1] 0 1cat("Number of fly genes")## Number of fly genesprint(dim(dtruth)[1])## [1] 13937cat("Number of edited genes")## Number of edited genesprint( dim(dtruth[ds_status==1,])[1] )## [1] 1000dtg <- dtruth[ds_status==1, gene] # list of tweaked genes in fly
dcg <- dtruth[ds_status==0, gene] # list of control genes in fly
cat("Values of Status")## Values of Statusprint( htruth[, unique(ds_status)] )## [1] 0 1cat("Number of human genes")## Number of human genesprint(dim(htruth)[1])## [1] 20410cat("Number of edited genes")## Number of edited genesprint( dim(htruth[ds_status==1,])[1] )## [1] 1000htg <- htruth[ds_status==1, gene] # list of tweaked genes in human
hcg <- htruth[ds_status==0, gene] # list of control genes in human
cat("Number of genes not edited in fly and human respectively")## Number of genes not edited in fly and human respectivelydxg <- rdk$Genes[!parent_id %in% dtruth$gene, parent_id]
print(length(dxg))## [1] 1745hxg <- rhk$Genes[!parent_id %in% htruth$gene, parent_id]
print(length(hxg))## [1] 41483These lists of edited and non-edited genes will serve as the basis for the evaluation of the tools. The simulated sets are restricted to protein-coding genes, whereas the quantifications were called on the full transcriptome annotations. Therefore, one would expect the expression levels of the excluded non-coding genes to be 0.
Not much info is given about the expression level and effect size of the genes making up the truth. We can use the abundances stored in the RATs output to get an idea of how the truth is distributed.
Flag truth status in DTU outputs
# Add flag fields, to mark the rows that are genes among the selected 1000 or the excluded.
# This prevents repeated subsetting and lookup later on.
# DRIMSeq
ddk[, selected := gene_id %in% dtg]
dds[, selected := gene_id %in% dtg]
dhk[, selected := gene_id %in% htg]
dhs[, selected := gene_id %in% htg]
ddk[, excluded := gene_id %in% dxg]
dds[, excluded := gene_id %in% dxg]
dhk[, excluded := gene_id %in% hxg]
dhs[, excluded := gene_id %in% hxg]
ddkt[, selected := gene_id %in% dtg]
ddst[, selected := gene_id %in% dtg]
dhkt[, selected := gene_id %in% htg]
dhst[, selected := gene_id %in% htg]
ddkt[, excluded := gene_id %in% dxg]
ddst[, excluded := gene_id %in% dxg]
dhkt[, excluded := gene_id %in% hxg]
dhst[, excluded := gene_id %in% hxg]
# SUPPA2
sdk[, selected := gene_id %in% dtg]
sds[, selected := gene_id %in% dtg]
shk[, selected := gene_id %in% htg]
shs[, selected := gene_id %in% htg]
sdk[, excluded := gene_id %in% dxg]
sds[, excluded := gene_id %in% dxg]
shk[, excluded := gene_id %in% hxg]
shs[, excluded := gene_id %in% hxg]
# DEXSeq
xds[, selected := geneID %in% dtg]
xhs[, selected := geneID %in% htg]
xds[, excluded := geneID %in% dxg]
xhs[, excluded := geneID %in% hxg]
xsds[, selected := geneID %in% dtg]
xshs[, selected := geneID %in% htg]
xsds[, excluded := geneID %in% dxg]
xshs[, excluded := geneID %in% hxg]
# RATs
rdk$Genes[, selected := parent_id %in% dtg]
rds$Genes[, selected := parent_id %in% dtg]
rhk$Genes[, selected := parent_id %in% htg]
rhs$Genes[, selected := parent_id %in% htg]
rdk$Transcripts[, selected := parent_id %in% dtg]
rds$Transcripts[, selected := parent_id %in% dtg]
rhk$Transcripts[, selected := parent_id %in% htg]
rhs$Transcripts[, selected := parent_id %in% htg]
rdk$Genes[, excluded := parent_id %in% dxg]
rds$Genes[, excluded := parent_id %in% dxg]
rhk$Genes[, excluded := parent_id %in% hxg]
rhs$Genes[, excluded := parent_id %in% hxg]
rdk$Transcripts[, excluded := parent_id %in% dxg]
rds$Transcripts[, excluded := parent_id %in% dxg]
rhk$Transcripts[, excluded := parent_id %in% hxg]
rhs$Transcripts[, excluded := parent_id %in% hxg]
rdkth$Genes[, selected := parent_id %in% dtg]
rdsth$Genes[, selected := parent_id %in% dtg]
rhkth$Genes[, selected := parent_id %in% htg]
rhsth$Genes[, selected := parent_id %in% htg]
rdkth$Transcripts[, selected := parent_id %in% dtg]
rdsth$Transcripts[, selected := parent_id %in% dtg]
rhkth$Transcripts[, selected := parent_id %in% htg]
rhsth$Transcripts[, selected := parent_id %in% htg]
rdkth$Genes[, excluded := parent_id %in% dxg]
rdsth$Genes[, excluded := parent_id %in% dxg]
rhkth$Genes[, excluded := parent_id %in% hxg]
rhsth$Genes[, excluded := parent_id %in% hxg]
rdkth$Transcripts[, excluded := parent_id %in% dxg]
rdsth$Transcripts[, excluded := parent_id %in% dxg]
rhkth$Transcripts[, excluded := parent_id %in% hxg]
rhsth$Transcripts[, excluded := parent_id %in% hxg]
rdkth2$Genes[, selected := parent_id %in% dtg]
rdsth2$Genes[, selected := parent_id %in% dtg]
rhkth2$Genes[, selected := parent_id %in% htg]
rhsth2$Genes[, selected := parent_id %in% htg]
rdkth2$Transcripts[, selected := parent_id %in% dtg]
rdsth2$Transcripts[, selected := parent_id %in% dtg]
rhkth2$Transcripts[, selected := parent_id %in% htg]
rhsth2$Transcripts[, selected := parent_id %in% htg]
rdkth2$Genes[, excluded := parent_id %in% dxg]
rdsth2$Genes[, excluded := parent_id %in% dxg]
rhkth2$Genes[, excluded := parent_id %in% hxg]
rhsth2$Genes[, excluded := parent_id %in% hxg]
rdkth2$Transcripts[, excluded := parent_id %in% dxg]
rdsth2$Transcripts[, excluded := parent_id %in% dxg]
rhkth2$Transcripts[, excluded := parent_id %in% hxg]
rhsth2$Transcripts[, excluded := parent_id %in% hxg]# DRIMSeq
ddk10[, selected := gene_id %in% dtg]
dds10[, selected := gene_id %in% dtg]
dhk10[, selected := gene_id %in% htg]
dhs10[, selected := gene_id %in% htg]
ddk10[, excluded := gene_id %in% dxg]
dds10[, excluded := gene_id %in% dxg]
dhk10[, excluded := gene_id %in% hxg]
dhs10[, excluded := gene_id %in% hxg]
ddkt10[, selected := gene_id %in% dtg]
ddst10[, selected := gene_id %in% dtg]
dhkt10[, selected := gene_id %in% htg]
dhst10[, selected := gene_id %in% htg]
ddkt10[, excluded := gene_id %in% dxg]
ddst10[, excluded := gene_id %in% dxg]
dhkt10[, excluded := gene_id %in% hxg]
dhst10[, excluded := gene_id %in% hxg]
# SUPPA2
sdk10[, selected := gene_id %in% dtg]
sds10[, selected := gene_id %in% dtg]
shk10[, selected := gene_id %in% htg]
shs10[, selected := gene_id %in% htg]
sdk10[, excluded := gene_id %in% dxg]
sds10[, excluded := gene_id %in% dxg]
shk10[, excluded := gene_id %in% hxg]
shs10[, excluded := gene_id %in% hxg]
# DEXSeq
xds10[, selected := geneID %in% dtg]
xhs10[, selected := geneID %in% htg]
xds10[, excluded := geneID %in% dxg]
xhs10[, excluded := geneID %in% hxg]
xsds10[, selected := geneID %in% dtg]
xshs10[, selected := geneID %in% htg]
xsds10[, excluded := geneID %in% dxg]
xshs10[, excluded := geneID %in% hxg]
# RATs
rdk10$Genes[, selected := parent_id %in% dtg]
rds10$Genes[, selected := parent_id %in% dtg]
rhk10$Genes[, selected := parent_id %in% htg]
rhs10$Genes[, selected := parent_id %in% htg]
rdk10$Transcripts[, selected := parent_id %in% dtg]
rds10$Transcripts[, selected := parent_id %in% dtg]
rhk10$Transcripts[, selected := parent_id %in% htg]
rhs10$Transcripts[, selected := parent_id %in% htg]
rdk10$Genes[, excluded := parent_id %in% dxg]
rds10$Genes[, excluded := parent_id %in% dxg]
rhk10$Genes[, excluded := parent_id %in% hxg]
rhs10$Genes[, excluded := parent_id %in% hxg]
rdk10$Transcripts[, excluded := parent_id %in% dxg]
rds10$Transcripts[, excluded := parent_id %in% dxg]
rhk10$Transcripts[, excluded := parent_id %in% hxg]
rhs10$Transcripts[, excluded := parent_id %in% hxg]
rdkth10$Genes[, selected := parent_id %in% dtg]
rdsth10$Genes[, selected := parent_id %in% dtg]
rhkth10$Genes[, selected := parent_id %in% htg]
rhsth10$Genes[, selected := parent_id %in% htg]
rdkth10$Transcripts[, selected := parent_id %in% dtg]
rdsth10$Transcripts[, selected := parent_id %in% dtg]
rhkth10$Transcripts[, selected := parent_id %in% htg]
rhsth10$Transcripts[, selected := parent_id %in% htg]
rdkth10$Genes[, excluded := parent_id %in% dxg]
rdsth10$Genes[, excluded := parent_id %in% dxg]
rhkth10$Genes[, excluded := parent_id %in% hxg]
rhsth10$Genes[, excluded := parent_id %in% hxg]
rdkth10$Transcripts[, excluded := parent_id %in% dxg]
rdsth10$Transcripts[, excluded := parent_id %in% dxg]
rhkth10$Transcripts[, excluded := parent_id %in% hxg]
rhsth10$Transcripts[, excluded := parent_id %in% hxg]
rdkth210$Genes[, selected := parent_id %in% dtg]
rdsth210$Genes[, selected := parent_id %in% dtg]
rhkth210$Genes[, selected := parent_id %in% htg]
rhsth210$Genes[, selected := parent_id %in% htg]
rdkth210$Transcripts[, selected := parent_id %in% dtg]
rdsth210$Transcripts[, selected := parent_id %in% dtg]
rhkth210$Transcripts[, selected := parent_id %in% htg]
rhsth210$Transcripts[, selected := parent_id %in% htg]
rdkth210$Genes[, excluded := parent_id %in% dxg]
rdsth210$Genes[, excluded := parent_id %in% dxg]
rhkth210$Genes[, excluded := parent_id %in% hxg]
rhsth210$Genes[, excluded := parent_id %in% hxg]
rdkth210$Transcripts[, excluded := parent_id %in% dxg]
rdsth210$Transcripts[, excluded := parent_id %in% dxg]
rhkth210$Transcripts[, excluded := parent_id %in% hxg]
rhsth210$Transcripts[, excluded := parent_id %in% hxg]# DRIMSeq
ddk50[, selected := gene_id %in% dtg]
dds50[, selected := gene_id %in% dtg]
dhk50[, selected := gene_id %in% htg]
dhs50[, selected := gene_id %in% htg]
ddk50[, excluded := gene_id %in% dxg]
dds50[, excluded := gene_id %in% dxg]
dhk50[, excluded := gene_id %in% hxg]
dhs50[, excluded := gene_id %in% hxg]
ddkt50[, selected := gene_id %in% dtg]
ddst50[, selected := gene_id %in% dtg]
dhkt50[, selected := gene_id %in% htg]
dhst50[, selected := gene_id %in% htg]
ddkt50[, excluded := gene_id %in% dxg]
ddst50[, excluded := gene_id %in% dxg]
dhkt50[, excluded := gene_id %in% hxg]
dhst50[, excluded := gene_id %in% hxg]
# SUPPA2
sdk50[, selected := gene_id %in% dtg]
sds50[, selected := gene_id %in% dtg]
shk50[, selected := gene_id %in% htg]
shs50[, selected := gene_id %in% htg]
sdk50[, excluded := gene_id %in% dxg]
sds50[, excluded := gene_id %in% dxg]
shk50[, excluded := gene_id %in% hxg]
shs50[, excluded := gene_id %in% hxg]
# DEXSeq
xds50[, selected := geneID %in% dtg]
xhs50[, selected := geneID %in% htg]
xds50[, excluded := geneID %in% dxg]
xhs50[, excluded := geneID %in% hxg]
xsds50[, selected := geneID %in% dtg]
xshs50[, selected := geneID %in% htg]
xsds50[, excluded := geneID %in% dxg]
xshs50[, excluded := geneID %in% hxg]
# RATs
rdk50$Genes[, selected := parent_id %in% dtg]
rds50$Genes[, selected := parent_id %in% dtg]
rhk50$Genes[, selected := parent_id %in% htg]
rhs50$Genes[, selected := parent_id %in% htg]
rdk50$Transcripts[, selected := parent_id %in% dtg]
rds50$Transcripts[, selected := parent_id %in% dtg]
rhk50$Transcripts[, selected := parent_id %in% htg]
rhs50$Transcripts[, selected := parent_id %in% htg]
rdk50$Genes[, excluded := parent_id %in% dxg]
rds50$Genes[, excluded := parent_id %in% dxg]
rhk50$Genes[, excluded := parent_id %in% hxg]
rhs50$Genes[, excluded := parent_id %in% hxg]
rdk50$Transcripts[, excluded := parent_id %in% dxg]
rds50$Transcripts[, excluded := parent_id %in% dxg]
rhk50$Transcripts[, excluded := parent_id %in% hxg]
rhs50$Transcripts[, excluded := parent_id %in% hxg]
rdkth50$Genes[, selected := parent_id %in% dtg]
rdsth50$Genes[, selected := parent_id %in% dtg]
rhkth50$Genes[, selected := parent_id %in% htg]
rhsth50$Genes[, selected := parent_id %in% htg]
rdkth50$Transcripts[, selected := parent_id %in% dtg]
rdsth50$Transcripts[, selected := parent_id %in% dtg]
rhkth50$Transcripts[, selected := parent_id %in% htg]
rhsth50$Transcripts[, selected := parent_id %in% htg]
rdkth50$Genes[, excluded := parent_id %in% dxg]
rdsth50$Genes[, excluded := parent_id %in% dxg]
rhkth50$Genes[, excluded := parent_id %in% hxg]
rhsth50$Genes[, excluded := parent_id %in% hxg]
rdkth50$Transcripts[, excluded := parent_id %in% dxg]
rdsth50$Transcripts[, excluded := parent_id %in% dxg]
rhkth50$Transcripts[, excluded := parent_id %in% hxg]
rhsth50$Transcripts[, excluded := parent_id %in% hxg]
rdkth250$Genes[, selected := parent_id %in% dtg]
rdsth250$Genes[, selected := parent_id %in% dtg]
rhkth250$Genes[, selected := parent_id %in% htg]
rhsth250$Genes[, selected := parent_id %in% htg]
rdkth250$Transcripts[, selected := parent_id %in% dtg]
rdsth250$Transcripts[, selected := parent_id %in% dtg]
rhkth250$Transcripts[, selected := parent_id %in% htg]
rhsth250$Transcripts[, selected := parent_id %in% htg]
rdkth250$Genes[, excluded := parent_id %in% dxg]
rdsth250$Genes[, excluded := parent_id %in% dxg]
rhkth250$Genes[, excluded := parent_id %in% hxg]
rhsth250$Genes[, excluded := parent_id %in% hxg]
rdkth250$Transcripts[, excluded := parent_id %in% dxg]
rdsth250$Transcripts[, excluded := parent_id %in% dxg]
rhkth250$Transcripts[, excluded := parent_id %in% hxg]
rhsth250$Transcripts[, excluded := parent_id %in% hxg]Distribution of gene expression
Plot of aggregated expression across all isoforms for each gene. Categorised by coding genes selected for DTU creation, coding genes not selected, and non-coding genes that were excluded from the expression simulation.
# Gene expression level of genes selected to be true DTU events.
df <- data.table(expression= c(rdk$Transcripts[, unique(totalA), by=parent_id]$V1,
rds$Transcripts[, unique(totalA), by=parent_id]$V1,
rhk$Transcripts[, unique(totalA), by=parent_id]$V1,
rhs$Transcripts[, unique(totalA), by=parent_id]$V1),
selected= c(rdk$Transcripts[, unique(selected), by=parent_id]$V1,
rds$Transcripts[, unique(selected), by=parent_id]$V1,
rhk$Transcripts[, unique(selected), by=parent_id]$V1,
rhs$Transcripts[, unique(selected), by=parent_id]$V1),
excluded= c(rdk$Transcripts[, unique(excluded), by=parent_id]$V1,
rds$Transcripts[, unique(excluded), by=parent_id]$V1,
rhk$Transcripts[, unique(excluded), by=parent_id]$V1,
rhs$Transcripts[, unique(excluded), by=parent_id]$V1),
cond= c(rep('Fruitfly', dim(rdk$Genes)[1] + dim(rds$Genes)[1]),
rep('Human', dim(rhk$Genes)[1] + dim(rhs$Genes)[1])),
quant= c(rep('Kallisto', dim(rdk$Genes)[1]), rep('Salmon', dim(rds$Genes)[1]),
rep('Kallisto', dim(rhk$Genes)[1]), rep('Salmon', dim(rhs$Genes)[1])))
df$SelectedCoding <- paste(as.character(df$selected), as.character(!df$excluded))
p <- ggplot(df, aes(x=expression, fill= SelectedCoding)) + gth +
facet_grid(cond ~ quant, scales='fixed') +
geom_histogram(position=position_dodge()) +
scale_fill_manual(values=c('gold', 'firebrick1', 'steelblue')) +
labs(x='total # isoform-length normalised reads', y='gene count', title='Gene expression levels in truth subset')
print(p)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#print(p + scale_x_continuous(limits=c(NA,100000)) + coord_cartesian(ylim=c(0,2000)))
print(p + scale_x_continuous(limits=c(NA,15000)) + coord_cartesian(ylim=c(0,500)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 4020 rows containing non-finite values (stat_bin).## Warning: Removed 8 rows containing missing values (geom_bar).
print(p + scale_x_continuous(limits=c(NA,2500)) + coord_cartesian(ylim=c(0,20)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 21107 rows containing non-finite values (stat_bin).
## Warning: Removed 8 rows containing missing values (geom_bar).
- “High expression” translates into mainly genes with at least ~500 reads.
- The shape of the distribution is very similar between the quantification tools and datasets, but the reasults from the two tools are not identical at higher magnification.
- The high expression genes may still contain non-expressed isoforms, so plotting the distribution of isoform expression without knowing which specific isoforms were selected within each selected gene would be uninformative.
- Both tools misassign reads to non-coding genes, from which no reads were simulated. The error is more pronounced for the human dataset, where there are many more non-coding annotated genes. Kallisto seems to be a little worse in this than Salmon for the human dataset, but does a little better than Salmon for the fruitfly dataset.
Distribution of effect size in the truth
How’s the effect size distributed between the genes selected for abundance editing, versus all the remaining genes (including the genes excluded from the simulation of abundances).
# Obtain FX size from RATs output.
# Since pre-filtering does not physically remove entries from RATs
# we have to use the `elig` flag to plot the data the meets the threshold.
df <- rbind(data.table(expression= c(rdk$Transcripts[(elig), Dprop],
rds$Transcripts[(elig), Dprop],
rhk$Transcripts[(elig), Dprop],
rhs$Transcripts[(elig), Dprop]),
selected= c(rdk$Transcripts[(elig), selected],
rds$Transcripts[(elig), selected],
rhk$Transcripts[(elig), selected],
rhs$Transcripts[(elig), selected]),
excluded= c(rdk$Transcripts[(elig), excluded],
rds$Transcripts[(elig), excluded],
rhk$Transcripts[(elig), excluded],
rhs$Transcripts[(elig), excluded]),
cond= c(rep('Fruitfly', dim(rdk$Transcripts[(elig),])[1] + dim(rds$Transcripts[(elig),])[1]),
rep('Human', dim(rhk$Transcripts[(elig),])[1] + dim(rhs$Transcripts[(elig),])[1])),
quant= c(rep('Kallisto', dim(rdk$Transcripts[(elig),])[1]), rep('Salmon', dim(rds$Transcripts[(elig),])[1]),
rep('Kallisto', dim(rhk$Transcripts[(elig),])[1]), rep('Salmon', dim(rhs$Transcripts[(elig),])[1])),
prefilter=0),
data.table(expression= c(rdk10$Transcripts[(elig), Dprop],
rds10$Transcripts[(elig), Dprop],
rhk10$Transcripts[(elig), Dprop],
rhs10$Transcripts[(elig), Dprop]),
selected= c(rdk10$Transcripts[(elig), selected],
rds10$Transcripts[(elig), selected],
rhk10$Transcripts[(elig), selected],
rhs10$Transcripts[(elig), selected]),
excluded= c(rdk10$Transcripts[(elig), excluded],
rds10$Transcripts[(elig), excluded],
rhk10$Transcripts[(elig), excluded],
rhs10$Transcripts[(elig), excluded]),
cond= c(rep('Fruitfly', dim(rdk10$Transcripts[(elig),])[1] + dim(rds10$Transcripts[(elig),])[1]),
rep('Human', dim(rhk10$Transcripts[(elig),])[1] + dim(rhs10$Transcripts[(elig),])[1])),
quant= c(rep('Kallisto', dim(rdk10$Transcripts[(elig),])[1]), rep('Salmon', dim(rds10$Transcripts[(elig),])[1]),
rep('Kallisto', dim(rhk10$Transcripts[(elig),])[1]), rep('Salmon', dim(rhs10$Transcripts[(elig),])[1])),
prefilter=10),
data.table(expression= c(rdk50$Transcripts[(elig), Dprop],
rds50$Transcripts[(elig), Dprop],
rhk50$Transcripts[(elig), Dprop],
rhs50$Transcripts[(elig), Dprop]),
selected= c(rdk50$Transcripts[(elig), selected],
rds50$Transcripts[(elig), selected],
rhk50$Transcripts[(elig), selected],
rhs50$Transcripts[(elig), selected]),
excluded= c(rdk50$Transcripts[(elig), excluded],
rds50$Transcripts[(elig), excluded],
rhk50$Transcripts[(elig), excluded],
rhs50$Transcripts[(elig), excluded]),
cond= c(rep('Fruitfly', dim(rdk50$Transcripts[(elig),])[1] + dim(rds50$Transcripts[(elig),])[1]),
rep('Human', dim(rhk50$Transcripts[(elig),])[1] + dim(rhs50$Transcripts[(elig),])[1])),
quant= c(rep('Kallisto', dim(rdk50$Transcripts[(elig),])[1]), rep('Salmon', dim(rds50$Transcripts[(elig),])[1]),
rep('Kallisto', dim(rhk50$Transcripts[(elig),])[1]), rep('Salmon', dim(rhs50$Transcripts[(elig),])[1])),
prefilter=50)
)
df$SelectedCoding <- paste(as.character(df$selected), as.character(!df$excluded))
p <- ggplot(df, aes(x = expression, fill = SelectedCoding)) + gth +
facet_grid(cond ~ quant, scales='fixed') +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_histogram(position = position_dodge()) +
scale_fill_manual(values=c('FALSE FALSE'='gold', 'FALSE TRUE'='firebrick1', 'TRUE TRUE'='steelblue')) +
labs(x='Dprop', title='RATs effect size distribution', subtitle="Abundance prefilter: {closest_state}")
anim_save("analysis_sim_v3_truthFXdist_rats.gif", p)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.anim_save("analysis_sim_v3_truthFXdist_rats1.gif", p + coord_cartesian(ylim=c(0,10000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.anim_save("analysis_sim_v3_truthFXdist_rats2.gif", p + coord_cartesian(ylim=c(0,1000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Are there really impossible Dprop values below -1? (there are NA Dprop values, so the outcome of this query will be either TRUE or NA)
print( any(rhs$Transcripts$Dprop < -1) )## [1] NA# Obtain FX size from SUPPA2 output.
df <- rbind(data.table(expression= c(sdk[, dpsi], sds[, dpsi], shk[, dpsi], shs[, dpsi]),
selected= c(sdk[, selected], sds[, selected], shk[, selected], shs[, selected]),
excluded= c(sdk[, excluded], sds[, excluded], shk[, excluded], shs[, excluded]),
cond= c(rep('Fruitfly', dim(sdk)[1] + dim(sds)[1]), rep('Human', dim(shk)[1] + dim(shs)[1])),
quant= c(rep('Kallisto', dim(sdk)[1]), rep('Salmon', dim(sds)[1]), rep('Kallisto', dim(shk)[1]), rep('Salmon', dim(shs)[1])),
prefilter=0),
data.table(expression= c(sdk10[, dpsi], sds10[, dpsi], shk10[, dpsi], shs10[, dpsi]),
selected= c(sdk10[, selected], sds10[, selected], shk10[, selected], shs10[, selected]),
excluded= c(sdk10[, excluded], sds10[, excluded], shk10[, excluded], shs10[, excluded]),
cond= c(rep('Fruitfly', dim(sdk10)[1] + dim(sds10)[1]), rep('Human', dim(shk10)[1] + dim(shs10)[1])),
quant= c(rep('Kallisto', dim(sdk10)[1]), rep('Salmon', dim(sds10)[1]), rep('Kallisto', dim(shk10)[1]), rep('Salmon', dim(shs10)[1])),
prefilter=10),
data.table(expression= c(sdk50[, dpsi], sds50[, dpsi], shk50[, dpsi], shs50[, dpsi]),
selected= c(sdk50[, selected], sds50[, selected], shk50[, selected], shs50[, selected]),
excluded= c(sdk50[, excluded], sds50[, excluded], shk50[, excluded], shs50[, excluded]),
cond= c(rep('Fruitfly', dim(sdk50)[1] + dim(sds50)[1]), rep('Human', dim(shk50)[1] + dim(shs50)[1])),
quant= c(rep('Kallisto', dim(sdk50)[1]), rep('Salmon', dim(sds50)[1]), rep('Kallisto', dim(shk50)[1]), rep('Salmon', dim(shs50)[1])),
prefilter=50)
)
df$SelectedCoding <- paste(as.character(df$selected), as.character(!df$excluded))
p <- ggplot(df, aes(x = expression, fill = SelectedCoding)) + gth +
facet_grid(cond ~ quant, scales='fixed') +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_histogram(position=position_dodge()) +
scale_fill_manual(values=c('FALSE FALSE'='gold', 'FALSE TRUE'='firebrick1', 'TRUE TRUE'='steelblue')) +
labs(x='dPSI', title='SUPPA2 effect size distribution', subtitle="Abundance prefilter: {closest_state}")
anim_save("analysis_sim_v3_truthFXdist_suppa.gif", p)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 519294 rows containing non-finite values (stat_bin).anim_save("analysis_sim_v3_truthFXdist_suppa1.gif", p + coord_cartesian(ylim=c(0,10000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 519294 rows containing non-finite values (stat_bin).anim_save("analysis_sim_v3_truthFXdist_suppa2.gif", p + coord_cartesian(ylim=c(0,1000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 519294 rows containing non-finite values (stat_bin).
# Obtained FX size from DRIMSeq pre-comparison object.
df <- rbind(data.table(expression= c(ddkt[, Dprop], ddst[, Dprop], dhkt[, Dprop], dhst[, Dprop]),
selected= c(ddkt[, selected], ddst[, selected], dhkt[, selected],dhst[, selected]),
excluded= c(ddkt[, excluded], ddst[, excluded], dhkt[, excluded],dhst[, excluded]),
cond= c(rep('Fruitfly', dim(ddkt)[1] + dim(ddst)[1]), rep('Human', dim(dhkt)[1] + dim(dhst)[1])),
quant= c(rep('Kallisto', dim(ddkt)[1]), rep('Salmon', dim(ddst)[1]), rep('Kallisto', dim(dhkt)[1]), rep('Salmon', dim(dhst)[1])),
prefilter=0),
data.table(expression= c(ddkt10[, Dprop], ddst10[, Dprop], dhkt10[, Dprop], dhst10[, Dprop]),
selected= c(ddkt10[, selected], ddst10[, selected], dhkt10[, selected], dhst10[, selected]),
excluded= c(ddkt10[, excluded], ddst10[, excluded], dhkt10[, excluded], dhst10[, excluded]),
cond= c(rep('Fruitfly', dim(ddkt10)[1] + dim(ddst10)[1]), rep('Human', dim(dhkt10)[1] + dim(dhst10)[1])),
quant= c(rep('Kallisto', dim(ddkt10)[1]), rep('Salmon', dim(ddst10)[1]), rep('Kallisto', dim(dhkt10)[1]), rep('Salmon', dim(dhst10)[1])),
prefilter=10),
data.table(expression= c(ddkt50[, Dprop], ddst50[, Dprop], dhkt50[, Dprop], dhst50[, Dprop]),
selected= c(ddkt50[, selected], ddst50[, selected], dhkt50[, selected], dhst50[, selected]),
excluded= c(ddkt50[, excluded], ddst50[, excluded], dhkt50[, excluded], dhst50[, excluded]),
cond= c(rep('Fruitfly', dim(ddkt50)[1] + dim(ddst50)[1]), rep('Human', dim(dhkt50)[1] + dim(dhst50)[1])),
quant= c(rep('Kallisto', dim(ddkt50)[1]), rep('Salmon', dim(ddst50)[1]), rep('Kallisto', dim(dhkt50)[1]), rep('Salmon', dim(dhst50)[1])),
prefilter=50)
)
df$SelectedCoding <- paste(as.character(df$selected), as.character(!df$excluded))
p <- ggplot(df, aes(x = expression, fill = selected)) + gth +
facet_grid(cond ~ quant, scales='fixed') +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_histogram(position=position_dodge()) +
scale_fill_manual(values=c('FALSE FALSE'='gold', 'FALSE TRUE'='firebrick1', 'TRUE TRUE'='steelblue')) +
labs(x='Dprop', title='DRIMSeq effect size distribution', subtitle="Abundance prefilter: {closest_state}")
anim_save("analysis_sim_v3_truthFXdist_drimseq.gif", p)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 666 rows containing non-finite values (stat_bin).anim_save("analysis_sim_v3_truthFXdist_drimseq1.gif", p + coord_cartesian(ylim=c(0,10000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 666 rows containing non-finite values (stat_bin).anim_save("analysis_sim_v3_truthFXdist_drimseq2.gif", p + coord_cartesian(ylim=c(0,1000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 666 rows containing non-finite values (stat_bin).
# Obtain FX size from DEXSeq output.
df <- rbind(data.table(expression= c(xds[, Dprop], xhs[, Dprop]),
selected= c(xds[, selected], xhs[, selected]),
excluded= c(xds[, excluded], xhs[, excluded]),
cond= c(rep('Fruitfly', dim(xds)[1]), rep('Human', dim(xhs)[1])),
quant= c(rep('Salmon', dim(xds)[1]), rep('Salmon', dim(xhs)[1])),
prefilter=0),
data.table(expression= c(xds10[, Dprop], xhs10[, Dprop]),
selected= c(xds10[, selected], xhs10[, selected]),
excluded= c(xds10[, excluded], xhs10[, excluded]),
cond= c(rep('Fruitfly', dim(xds10)[1]), rep('Human', dim(xhs10)[1])),
quant= c(rep('Salmon', dim(xds10)[1]), rep('Salmon', dim(xhs10)[1])),
prefilter=10),
data.table(expression= c(xds50[, Dprop], xhs50[, Dprop]),
selected= c(xds50[, selected], xhs50[, selected]),
excluded= c(xds10[, excluded], xhs10[, excluded]),
cond= c(rep('Fruitfly', dim(xds50)[1]), rep('Human', dim(xhs50)[1])),
quant= c(rep('Salmon', dim(xds50)[1]), rep('Salmon', dim(xhs50)[1])),
prefilter=50)
)
df$SelectedCoding <- paste(as.character(df$selected), as.character(!df$excluded))
p <- ggplot(df, aes(x = expression, fill = selected)) + gth +
facet_grid(cond ~ quant, scales='fixed') +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_histogram(position=position_dodge()) +
scale_fill_manual(values=c('FALSE FALSE'='gold', 'FALSE TRUE'='firebrick1', 'TRUE TRUE'='steelblue')) +
labs(x='Dprop', title='DEXSeq effect size distribution', subtitle="Abundance prefilter: {closest_state}")
anim_save("analysis_sim_v3_truthFXdist_dexseq.gif", p)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 290704 rows containing non-finite values (stat_bin).anim_save("analysis_sim_v3_truthFXdist_dexseq1.gif", p + coord_cartesian(ylim=c(0,10000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 290704 rows containing non-finite values (stat_bin).anim_save("analysis_sim_v3_truthFXdist_dexseq2.gif", p + coord_cartesian(ylim=c(0,1000)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 290704 rows containing non-finite values (stat_bin).
- The apparent existence of Dprop/dPSI values below -1, is an artifact of the binning method.
- The selected genes often contain more than 2 isoforms, but only 2 were edited to create DTU, so the high abundance of 0 effect size in selected genes is not surprising.
- There is a healthily wide distribution of non-zero effect sizes in the selected genes.
- The Dprop/dPSI distribution should be symmetrical around 0, by nature of how the datasets were created. Deviation from symmetry reflects the stochasticity in read simulation and quantification.
- Although for all four tools the input abundances are the same, there are small differences in the effect size distributions, possibly as a result of internal filtering, normalisation, or method of calculation.
- The effect sizes of negative genes (both coding and non-coding) appear wider distributed than one might hope, but even a 10-read pre-filter greatly reduces the number of high value occurencies.
- Higher abundance pre-filter values mainly reduce the number of negative genes (both coding and non-coding). However this is not informative for real-life datasets, as the positive genes were explicitly selected to have high expression, well above the maximum threshold used here as shown in the previous section.
Individual comparisons to truth
The truth consists of swapping abundance for the two most abundant isoforms of 1000 well-expressed genes in each dataset. Ideally, tools are expected to identify all of these 1000 genes or 2000 transcripts, and only those. As all genes were selected to have high expression, false negatives are expected only for very small effect sizes. No a-priori limit was used in generating the dataset, and no such pre-filter was applied prior to DTU calling. As all other genes and transcripts should be completely unchanged by design, all false positives must be attributed to random fluctuations in the simulation and quantification processes: The datasets were simulated to identical final abundances, except for the 1000 selected genes, but the reads were created independently for each set. Thus, the exact set of reads would influence the quantification.
We track the following values:
- TP - true positives (designed to be DTU and found to be DTU).
- FP - false positives (designed to be DTU but not found to be DTU).
- TN - true negatives (designed not to be DTU and found not to be DTU).
- FN - false negatives (designed not to be DTU but found to be DTU).
- TNA - an excluded (non-coding) gene (and therefore not “expressed”) that is found not to be DTU. This is a type of True Negative, in that we should expect no DTU for the genes excluded from the simulation.
- FNA - an excluded (non-coding) gene falsely found to be DTU. This is a type of False Positive, in that we should expect no DTU for the genes excluded from the simulation.
# Colour code
col <- c('TP'='green', 'FP'='red', 'TN'='blue', 'FN'='darkorange', 'TNA'='grey50', 'FNA'='pink')Comparison between the tools at given thresholds, ignoring genes excluded from the simulations.
SUPPA2 and DRIMSeq do not apply thresholds or classify DTU events. They return the test results raw. RATs also reports the raw results, but the reproducibility feature requires that specific significance, effect size and noise thresholds be defined in advance, and any post-hoc change of these thresholds invalidates the reproducibility measurements. However, the reproducibility thresholds can be adjusted post-hoc, with no need for a rerun.
The simulated sets were focused on coding genes only. Non-coding genes would be expected to appear as non-expressed and could be added to the FP and TN. Here, we opted to calculate the performance metrics using only the genes explicitly designated as either positive or negative.
Assistant functions
# Count TP, FP, TN, FN
# Not used directly.
countup <- function(x, tool, fx, qt=0.95, rt=0.85) {
if (tool=='RATs(g)') { # RATs Gene level
# Re-classify DTU
dk <- x[[1]]$Gene[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt]
ds <- x[[2]]$Gene[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt]
hk <- x[[3]]$Gene[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt]
hs <- x[[4]]$Gene[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt]
# Structure.
return( data.table(tp= c(dim( x[[1]]$Genes[dk & selected, TRUE])[1],
dim(x[[2]]$Genes[ds & selected, TRUE])[1],
dim(x[[3]]$Genes[hk & selected, TRUE])[1],
dim(x[[4]]$Genes[hs & selected, TRUE])[1] ),
fp= c(dim(x[[1]]$Genes[dk & !(selected | excluded), TRUE])[1],
dim(x[[2]]$Genes[ds & !(selected | excluded), TRUE])[1],
dim(x[[3]]$Genes[hk & !(selected | excluded), TRUE])[1],
dim(x[[4]]$Genes[hs & !(selected | excluded), TRUE])[1] ),
tn= c(dim(x[[1]]$Genes[(!dk) & !(selected | excluded), TRUE])[1],
dim(x[[2]]$Genes[(!ds) & !(selected | excluded), TRUE])[1],
dim(x[[3]]$Genes[(!hk) & !(selected | excluded), TRUE])[1],
dim(x[[4]]$Genes[(!hs) & !(selected | excluded), TRUE])[1] ),
fn= c(dim(x[[1]]$Genes[(!dk) & selected, TRUE])[1],
dim(x[[2]]$Genes[(!ds) & selected, TRUE])[1],
dim(x[[3]]$Genes[(!hk) & selected, TRUE])[1],
dim(x[[4]]$Genes[(!hs) & selected, TRUE])[1] ),
fna= c(dim(x[[1]]$Genes[dk & excluded, TRUE])[1],
dim(x[[2]]$Genes[ds & excluded, TRUE])[1],
dim(x[[3]]$Genes[hk & excluded, TRUE])[1],
dim(x[[4]]$Genes[hs & excluded, TRUE])[1] ),
tna= c(dim(x[[1]]$Genes[(!dk) & !excluded, TRUE])[1],
dim(x[[2]]$Genes[(!ds) & !excluded, TRUE])[1],
dim(x[[3]]$Genes[(!hk) & !excluded, TRUE])[1],
dim(x[[4]]$Genes[(!hs) & !excluded, TRUE])[1] ),
set= c('Fruitfly', 'Fruitfly', 'Human', 'Human'),
quant= c('Kallisto', 'Salmon', 'Kallisto', 'Salmon') ) )
} else if (tool=='RATs(t)') { # RATs transcript level
# Reclassify DTU.
tmp <- data.table(gene = x[[1]]$Transcript$parent_id, sel = x[[1]]$Transcript$selected, xcl = x[[1]]$Transcript$excluded,
thit = x[[1]]$Transcripts[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt], ghit = NA)
tmp[, ghit := any(thit), by= gene]
dk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[2]]$Transcript$parent_id, sel = x[[2]]$Transcript$selected, xcl = x[[2]]$Transcript$excluded,
thit = x[[2]]$Transcripts[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt], ghit = NA)
tmp[, ghit := any(thit), by= gene]
ds <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[3]]$Transcript$parent_id, sel = x[[3]]$Transcript$selected, xcl = x[[3]]$Transcript$excluded,
thit = x[[3]]$Transcript[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt], ghit = NA)
tmp[, ghit := any(thit), by= gene]
hk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[4]]$Transcript$parent_id, sel = x[[4]]$Transcript$selected, xcl = x[[4]]$Transcript$excluded,
thit = x[[4]]$Transcript[, !is.na(quant_dtu_freq) & quant_dtu_freq >= qt & !is.na(rep_dtu_freq) & rep_dtu_freq >= rt], ghit = NA)
tmp[, ghit := any(thit), by= gene]
hs <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
# Structure.
return( data.table(tp= c(dk['tp'], ds['tp'], hk['tp'], hs['tp']),
fp= c(dk['fp'], ds['fp'], hk['fp'], hs['fp']),
fn= c(dk['fn'], ds['fn'], hk['fn'], hs['fn']),
tn= c(dk['tn'], ds['tn'], hk['tn'], hs['tn']),
fna= c(dk['fna'], ds['fna'], hk['fna'], hs['fna']),
tna= c(dk['tna'], ds['tna'], hk['tna'], hs['tna']),
set= c('Fruitfly', 'Fruitfly', 'Human', 'Human'),
quant= c('Kallisto', 'Salmon', 'Kallisto', 'Salmon') ) )
} else if (tool=='SUPPA2') {
# Apply thresholds and aggregate by gene
tmp <- data.table(gene = x[[1]]$gene_id, sel = x[[1]]$selected, xcl = x[[1]]$excluded, thit = x[[1]]$dpsi >= fx & x[[1]]$pval < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
dk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[2]]$gene_id, sel = x[[2]]$selected, xcl = x[[2]]$excluded, thit = x[[2]]$dpsi >= fx & x[[2]]$pval < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
ds <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[3]]$gene_id, sel = x[[3]]$selected, xcl = x[[3]]$excluded, thit = x[[3]]$dpsi >= fx & x[[3]]$pval < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
hk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[4]]$gene_id, sel = x[[4]]$selected, xcl = x[[4]]$excluded, thit = x[[4]]$dpsi >= fx & x[[4]]$pval < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
hs <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
# Lack of report as FALSE, like RATs does.
dk[is.na(dk)] <- FALSE
ds[is.na(ds)] <- FALSE
hk[is.na(hk)] <- FALSE
hk[is.na(hs)] <- FALSE
# Structure.
return( data.table(tp= c(dk['tp'], ds['tp'], hk['tp'], hs['tp']),
fp= c(dk['fp'], ds['fp'], hk['fp'], hs['fp']),
fn= c(dk['fn'], ds['fn'], hk['fn'], hs['fn']),
tn= c(dk['tn'], ds['tn'], hk['tn'], hs['tn']),
fna= c(dk['fna'], ds['fna'], hk['fna'], hs['fna']),
tna= c(dk['tna'], ds['tna'], hk['tna'], hs['tna']),
set= c('Fruitfly', 'Fruitfly', 'Human', 'Human'),
quant= c('Kallisto', 'Salmon', 'Kallisto', 'Salmon') ) )
} else if (tool=='DRIMSeq(g)') {
# Apply thresholds. any() is used to reduce a group of identical values to just one.
dk <- x[[1]][, adj_pvalue < 0.05] & x[[1]][, maxDprop >= fx]
ds <- x[[2]][, adj_pvalue < 0.05] & x[[2]][, maxDprop >= fx]
hk <- x[[3]][, adj_pvalue < 0.05] & x[[3]][, maxDprop >= fx]
hs <- x[[4]][, adj_pvalue < 0.05] & x[[4]][, maxDprop >= fx]
# Lack of report as FALSE, like RATs does.
dk[is.na(dk)] <- FALSE
ds[is.na(ds)] <- FALSE
hk[is.na(hk)] <- FALSE
hk[is.na(hs)] <- FALSE
# Structure.
return( data.table(tp= c(dim(x[[1]][dk & selected, TRUE])[1],
dim(x[[2]][ds & selected, TRUE])[1],
dim(x[[3]][hk & selected, TRUE])[1],
dim(x[[4]][hs & selected, TRUE])[1] ),
fp= c(dim(x[[1]][dk & !(selected | excluded), TRUE])[1],
dim(x[[2]][ds & !(selected | excluded), TRUE])[1],
dim(x[[3]][hk & !(selected | excluded), TRUE])[1],
dim(x[[4]][hs & !(selected | excluded), TRUE])[1] ),
fn= c(dim(x[[1]][(!dk) & selected, TRUE])[1],
dim(x[[2]][(!ds) & selected, TRUE])[1],
dim(x[[3]][(!hk) & selected, TRUE])[1],
dim(x[[4]][(!hs) & selected, TRUE])[1] ),
tn= c(dim(x[[1]][(!dk) & !(selected | excluded), TRUE])[1],
dim(x[[2]][(!ds) & !(selected | excluded), TRUE])[1],
dim(x[[3]][(!hk) & !(selected | excluded), TRUE])[1],
dim(x[[4]][(!hs) & !(selected | excluded), TRUE])[1] ),
fna= c(dim(x[[1]][dk & excluded, TRUE])[1],
dim(x[[2]][ds & excluded, TRUE])[1],
dim(x[[3]][hk & excluded, TRUE])[1],
dim(x[[4]][hs & excluded, TRUE])[1] ),
tna= c(dim(x[[1]][(!dk) & !excluded, TRUE])[1],
dim(x[[2]][(!ds) & !excluded, TRUE])[1],
dim(x[[3]][(!hk) & !excluded, TRUE])[1],
dim(x[[4]][(!hs) & !excluded, TRUE])[1] ),
set= c('Fruitfly', 'Fruitfly', 'Human', 'Human'),
quant= c('Kallisto', 'Salmon', 'Kallisto', 'Salmon') ) )
} else if (tool=='DRIMSeq(t)') {
# Apply thresholds and aggregate by gene
tmp <- data.table(gene = x[[1]]$gene_id, sel = x[[1]]$selected, xcl = x[[1]]$excluded, thit = x[[1]]$Dprop >= fx & x[[1]]$adj_pvalue < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
dk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[2]]$gene_id, sel = x[[2]]$selected, xcl = x[[2]]$excluded, thit = x[[2]]$Dprop >= fx & x[[2]]$adj_pvalue < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
ds <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[3]]$gene_id, sel = x[[3]]$selected, xcl = x[[3]]$excluded, thit = x[[3]]$Dprop >= fx & x[[3]]$adj_pvalue < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
hk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[4]]$gene_id, sel = x[[4]]$selected, xcl = x[[4]]$excluded, thit = x[[4]]$Dprop >= fx & x[[4]]$adj_pvalue < 0.05, ghit = NA)
tmp[, ghit := any(thit), by= gene]
hs <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
# Lack of report as FALSE, like RATs does.
dk[is.na(dk)] <- FALSE
ds[is.na(ds)] <- FALSE
hk[is.na(hk)] <- FALSE
hk[is.na(hs)] <- FALSE
# Structure.
return( data.table(tp= c(dk['tp'], ds['tp'], hk['tp'], hs['tp']),
fp= c(dk['fp'], ds['fp'], hk['fp'], hs['fp']),
fn= c(dk['fn'], ds['fn'], hk['fn'], hs['fn']),
tn= c(dk['tn'], ds['tn'], hk['tn'], hs['tn']),
fna= c(dk['fna'], ds['fna'], hk['fna'], hs['fna']),
tna= c(dk['tna'], ds['tna'], hk['tna'], hs['tna']),
set= c('Fruitfly', 'Fruitfly', 'Human', 'Human'),
quant= c('Kallisto', 'Salmon', 'Kallisto', 'Salmon') ) )
} else if (tool=='DEXSeq') {
# Apply thresholds and aggregate by gene
tmp <- data.table(gene = x[[1]]$geneID, sel = x[[1]]$selected, xcl = x[[1]]$excluded, thit = x[[1]]$padj < 0.05 & x[[1]]$Dprop >= fx, ghit = NA)
tmp[is.na(thit), thit := FALSE] # Ineligible genes are NA, but we'll count them as FALSE, just like RATs does.
tmp[, ghit := any(thit), by= gene]
dk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[2]]$geneID, sel = x[[2]]$selected, xcl = x[[2]]$excluded, thit = x[[2]]$transcript < 0.05 & x[[2]]$Dprop >= fx, ghit = NA)
tmp[is.na(thit), thit := FALSE]
tmp[, ghit := any(thit), by= gene]
ds <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[3]]$geneID, sel = x[[3]]$selected, xcl = x[[3]]$excluded, thit = x[[3]]$padj < 0.05 & x[[3]]$Dprop >= fx, ghit = NA)
tmp[is.na(thit), thit := FALSE]
tmp[, ghit := any(thit), by= gene]
hk <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
tmp <- data.table(gene = x[[4]]$geneID, sel = x[[4]]$selected, xcl = x[[4]]$excluded, thit = x[[4]]$transcript < 0.05 & x[[4]]$Dprop >= fx, ghit = NA)
tmp[is.na(thit), thit := FALSE]
tmp[, ghit := any(thit), by= gene]
hs <- c( tp = length(tmp[ghit & sel, TRUE, by=gene]$V1),
fp = length(tmp[ghit & !(sel | xcl), TRUE, by=gene]$V1),
fn = length(tmp[(!ghit) & sel, TRUE, by=gene]$V1),
tn = length(tmp[(!ghit) & !(sel | xcl), TRUE, by=gene]$V1),
fna = length(tmp[(ghit) & xcl, TRUE, by=gene]$V1),
tna = length(tmp[(!ghit) & xcl, TRUE, by=gene]$V1) )
# Lack of report as FALSE, like RATs does.
dk[is.na(dk)] <- FALSE
ds[is.na(ds)] <- FALSE
hk[is.na(hk)] <- FALSE
hk[is.na(hs)] <- FALSE
# Structure.
return( data.table(tp= c(dk['tp'], ds['tp'], hk['tp'], hs['tp']),
fp= c(dk['fp'], ds['fp'], hk['fp'], hs['fp']),
fn= c(dk['fn'], ds['fn'], hk['fn'], hs['fn']),
tn= c(dk['tn'], ds['tn'], hk['tn'], hs['tn']),
fna= c(dk['fna'], ds['fna'], hk['fna'], hs['fna']),
tna= c(dk['tna'], ds['tna'], hk['tna'], hs['tna']),
set= c('Fruitfly', 'Fruitfly', 'Human', 'Human'),
quant= c('Salmon', 'Salmon', 'Salmon', 'Salmon') ) )
} else {
stop('Tool?')
}
}
# Calculate performance metrics.
# x: List in order fly/kallisto, fly/salmon, human/kallisto, human/salmon
measureup <- function(x, tool, fx, qt=0.95, rt=0.85, pre=0) {
df <- countup(x, tool, fx, qt, rt)
dt <- data.table(value= c(round(df[, tp/(tp+fn)], digits=3),
round(df[, fp/(tp+fp)], digits=3),
round(df[, as.double(tp*tn - fp*fn) / sqrt( as.double(tp + fp)*as.double(tp + fn)*as.double(tn + fp)*as.double(tn + fn) )], digits=3)),
type= factor(c(rep('Sensitivity', 4), rep('FDR', 4), rep('MCC', 4)), levels=c('Sensitivity', 'FDR', 'MCC')),
tool= factor(rep(tool, 12), levels=c('DRIMSeq(g)', 'DRIMSeq(t)', 'RATs(g)', 'RATs(t)', 'SUPPA2', 'DEXSeq')),
species=rep(df$set, 3),
quantification=rep(df$quant, 3),
effect=rep(paste('D >= ', fx), 12) ,
reproducibility=factor(rep(paste('Qrep', ifelse(is.na(qt), '', '>='), qt, 'Rrep', ifelse(is.na(rt), '', '>='), rt), 12)),
prefilter=pre)
if (tool=='DEXSeq')
dt[, tool := rep(c('DEXSeq', 'DEXsgr'), 6)]
return(dt)
}Performance at pval < 0.05
RATs
# RATs Dprop >= 0.20
rgd <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.95, rt=0.85, pre=0)
rtd <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.95, rt=0.85, pre=0)
# more permissive reproducibility
rgd1 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.90, rt=0.75, pre=0)
rtd1 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.90, rt=0.75, pre=0)
rgd2 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.85, rt=0.65, pre=0)
rtd2 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.85, rt=0.65, pre=0)
rgd3 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.80, rt=0.55, pre=0)
rtd3 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.80, rt=0.55, pre=0)
# RATs Dprop >= 0.10
rgf <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.95, rt=0.85, pre=0)
rtf <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.95, rt=0.85, pre=0)
# more permissive reproducibility
rgf1 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.90, rt=0.75, pre=0)
rtf1 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.90, rt=0.75, pre=0)
rgf2 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.85, rt=0.65, pre=0)
rtf2 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.85, rt=0.65, pre=0)
rgf3 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.80, rt=0.55, pre=0)
rtf3 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.80, rt=0.55, pre=0)
# RATs Dprop >= 0.05
rgff <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.95, rt=0.85, pre=0)
rtff <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.95, rt=0.85, pre=0)
# more permissive reproducibility
rgff1 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.90, rt=0.75, pre=0)
rtff1 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.90, rt=0.75, pre=0)
rgff2 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.85, rt=0.65, pre=0)
rtff2 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.85, rt=0.65, pre=0)
rgff3 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.80, rt=0.55, pre=0)
rtff3 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.80, rt=0.55, pre=0)# RATs Dprop >= 0.20
rgd10 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.95, rt=0.85, pre=10)
rtd10 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.95, rt=0.85, pre=10)
# more permissive reproducibility
rgd110 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.90, rt=0.75, pre=10)
rtd110 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.90, rt=0.75, pre=10)
rgd210 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.85, rt=0.65, pre=10)
rtd210 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.85, rt=0.65, pre=10)
rgd310 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.80, rt=0.55, pre=10)
rtd310 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.80, rt=0.55, pre=10)
# RATs Dprop >= 0.10
rgf10 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.95, rt=0.85, pre=10)
rtf10 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.95, rt=0.85, pre=10)
# more permissive reproducibility
rgf110 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.90, rt=0.75, pre=10)
rtf110 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.90, rt=0.75, pre=10)
rgf210 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.85, rt=0.65, pre=10)
rtf210 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.85, rt=0.65, pre=10)
rgf310 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.80, rt=0.55, pre=10)
rtf310 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.80, rt=0.55, pre=10)
# RATs Dprop >= 0.05
rgff10 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.95, rt=0.85, pre=10)
rtff10 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.95, rt=0.85, pre=10)
# more permissive reproducibility
rgff110 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.90, rt=0.75, pre=10)
rtff110 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.90, rt=0.75, pre=10)
rgff210 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.85, rt=0.65, pre=10)
rtff210 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.85, rt=0.65, pre=10)
rgff310 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.80, rt=0.55, pre=10)
rtff310 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.80, rt=0.55, pre=10)# RATs Dprop >= 0.20
rgd50 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.95, rt=0.85, pre=50)
rtd50 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.95, rt=0.85, pre=50)
# more permissive reproducibility
rgd150 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.90, rt=0.75, pre=50)
rtd150 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.90, rt=0.75, pre=50)
rgd250 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.85, rt=0.65, pre=50)
rtd250 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.85, rt=0.65, pre=50)
rgd350 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.80, rt=0.55, pre=50)
rtd350 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.80, rt=0.55, pre=50)
# RATs Dprop >= 0.10
rgf50 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.95, rt=0.85, pre=50)
rtf50 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.95, rt=0.85, pre=50)
# more permissive reproducibility
rgf150 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.90, rt=0.75, pre=50)
rtf150 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.90, rt=0.75, pre=50)
rgf250 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.85, rt=0.65, pre=50)
rtf250 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.85, rt=0.65, pre=50)
rgf350 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.80, rt=0.55, pre=50)
rtf350 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.80, rt=0.55, pre=50)
# RATs Dprop >= 0.05
rgff50 <- measureup(list(rdkth250,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.95, rt=0.85, pre=50)
rtff50 <- measureup(list(rdkth250,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.95, rt=0.85, pre=50)
# more permissive reproducibility
rgff150 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(g)', fx=0.05, qt=0.90, rt=0.75, pre=50)
rtff150 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(t)', fx=0.05, qt=0.90, rt=0.75, pre=50)
rgff250 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(g)', fx=0.05, qt=0.85, rt=0.65, pre=50)
rtff250 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(t)', fx=0.05, qt=0.85, rt=0.65, pre=50)
rgff350 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(g)', fx=0.05, qt=0.80, rt=0.55, pre=50)
rtff350 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(t)', fx=0.05, qt=0.80, rt=0.55, pre=50)SUPPA2
# SUPPA2 does not bootstrap the reproducibility.
# SUPPA2 dPSI >= 0.20
sd <- measureup(list(sdk,sds,shk,shs), 'SUPPA2', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0)
# SUPPA2 dPSI >= 0.10
sf <- measureup(list(sdk,sds,shk,shs), 'SUPPA2', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0)
# SUPPA2 dPSI >= 0.05
sff <- measureup(list(sdk,sds,shk,shs), 'SUPPA2', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0)# SUPPA2 does not bootstrap the reproducibility.
# SUPPA2 dPSI >= 0.20
sd10 <- measureup(list(sdk10,sds10,shk10,shs10), 'SUPPA2', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10)
# SUPPA2 dPSI >= 0.10
sf10 <- measureup(list(sdk10,sds10,shk10,shs10), 'SUPPA2', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10)
# SUPPA2 dPSI >= 0.05
sff10 <- measureup(list(sdk10,sds10,shk10,shs10), 'SUPPA2', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10)# SUPPA2 does not bootstrap the reproducibility.
# SUPPA2 dPSI >= 0.20
sd50 <- measureup(list(sdk50,sds50,shk50,shs50), 'SUPPA2', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50)
# SUPPA2 dPSI >= 0.10
sf50 <- measureup(list(sdk50,sds50,shk50,shs50), 'SUPPA2', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50)
# SUPPA2 dPSI >= 0.05
sff50 <- measureup(list(sdk50,sds50,shk50,shs50), 'SUPPA2', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50)DRIMSeq
# DRIMSeq does not bootstrap the reproducibility.
# Gene-level
# Dprop >= 0.20
dd <- measureup(list(ddk,dds,dhk,dhs), 'DRIMSeq(g)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.10
df <- measureup(list(ddk,dds,dhk,dhs), 'DRIMSeq(g)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.05
dff <- measureup(list(ddk,dds,dhk,dhs), 'DRIMSeq(g)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0)
# Transcript-level
# Dprop >= 0.20
ddt <- measureup(list(ddkt,ddst,dhkt,dhst), 'DRIMSeq(t)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.10
dft <- measureup(list(ddkt,ddst,dhkt,dhst), 'DRIMSeq(t)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.05
dfft <- measureup(list(ddkt,ddst,dhkt,dhst), 'DRIMSeq(t)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0)# DRIMSeq does not bootstrap the reproducibility.
# Gene-level
# Dprop >= 0.20
dd10 <- measureup(list(ddk10,dds10,dhk10,dhs10), 'DRIMSeq(g)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.10
df10 <- measureup(list(ddk10,dds10,dhk10,dhs10), 'DRIMSeq(g)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.05
dff10 <- measureup(list(ddk10,dds10,dhk10,dhs10), 'DRIMSeq(g)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10)
# Transcript-level
# Dprop >= 0.20
ddt10 <- measureup(list(ddkt10,ddst10,dhkt10,dhst10), 'DRIMSeq(t)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.10
dft10 <- measureup(list(ddkt10,ddst10,dhkt10,dhst10), 'DRIMSeq(t)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.05
dfft10 <- measureup(list(ddkt10,ddst10,dhkt10,dhst10), 'DRIMSeq(t)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10)# DRIMSeq does not bootstrap the reproducibility.
# Gene-level
# Dprop >= 0.20
dd50 <- measureup(list(ddk50,dds50,dhk50,dhs50), 'DRIMSeq(g)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.10
df50 <- measureup(list(ddk50,dds50,dhk50,dhs50), 'DRIMSeq(g)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.05
dff50 <- measureup(list(ddk50,dds50,dhk50,dhs50), 'DRIMSeq(g)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50)
# Transcript-level
# Dprop >= 0.20
ddt50 <- measureup(list(ddkt50,ddst50,dhkt50,dhst50), 'DRIMSeq(t)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.10
dft50 <- measureup(list(ddkt50,ddst50,dhkt50,dhst50), 'DRIMSeq(t)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.05
dfft50 <- measureup(list(ddkt50,ddst50,dhkt50,dhst50), 'DRIMSeq(t)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50)DEXSeq
# DEXSeq does not bootstrap the reproducibility.
# DEXSeq with and without stageR corrections. Dprop >= 0.20
xd <- measureup(list(xds, xsds, xhs, xshs), 'DEXSeq', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0 )
# DEXSeq with and without stageR corrections. Dprop >= 0.10
xf <- measureup(list(xds, xsds, xhs, xshs), 'DEXSeq', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0 )
# DEXSeq with and without stageR corrections. Dprop >= 0.05
xff <- measureup(list(xds, xsds, xhs, xshs), 'DEXSeq', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0 )# DEXSeq does not bootstrap the reproducibility.
# DEXSeq with and without stageR corrections. Dprop >= 0.20
xd10 <- measureup(list(xds10, xsds10, xhs10, xshs10), 'DEXSeq', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10 )
# DEXSeq with and without stageR corrections. Dprop >= 0.10
xf10 <- measureup(list(xds10, xsds10, xhs10, xshs10), 'DEXSeq', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10 )
# DEXSeq with and without stageR corrections. Dprop >= 0.05
xff10 <- measureup(list(xds10, xsds10, xhs10, xshs10), 'DEXSeq', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10 )# DEXSeq does not bootstrap the reproducibility.
# DEXSeq with and without stageR corrections. Dprop >= 0.20
xd50 <- measureup(list(xds50, xsds50, xhs50, xshs50), 'DEXSeq', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50 )
# DEXSeq with and without stageR corrections. Dprop >= 0.10
xf50 <- measureup(list(xds50, xsds50, xhs50, xshs50), 'DEXSeq', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50 )
# DEXSeq with and without stageR corrections. Dprop >= 0.05
xff50 <- measureup(list(xds50, xsds50, xhs50, xshs50), 'DEXSeq', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50 )Plot
Merge results into single table.
# Easier to work with ggplot2 if everything is in one table.
mega <- rbindlist(list(rgd, rtd, rgd1, rtd1, rgd2, rtd2, rgd3, rtd3,
rgf, rtf, rgf1, rtf1, rgf2, rtf2, rgf3, rtf3,
rgff, rtff, rgff1, rtff1, rgff2, rtff2, rgff3, rtff3,
sd, sf, sff, dd, df, dff, ddt, dft, dfft, xd, xf, xff,
rgd10, rtd10, rgd110, rtd110, rgd210, rtd210, rgd310, rtd310,
rgf10, rtf10, rgf110, rtf110, rgf210, rtf210, rgf310, rtf310,
rgff10, rtff10, rgff110, rtff110, rgff210, rtff210, rgff310, rtff310,
sd10, sf10, sff10, dd10, df10, dff10, ddt10, dft10, dfft10, xd10, xf10, xff10,
rgd50, rtd50, rgd150, rtd150, rgd250, rtd250, rgd350, rtd350,
rgf50, rtf50, rgf150, rtf150, rgf250, rtf250, rgf350, rtf350,
rgff50, rtff50, rgff150, rtff150, rgff250, rtff250, rgff350, rtff350,
sd50, sf50, sff50, dd50, df50, dff50, ddt50, dft50, dfft50, xd50, xf50, xff50))Plot.
# Reduce the number of points to include in the MS figure.
ms <- mega[! mega$reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65'), ]
ms <- ms[quantification == 'Salmon', ]
ms <- ms[prefilter == 0 | prefilter == 50]
ms <- ms[tool != 'DEXsgr'] # stageR is a post-processing that could be applied to any tool. For fairness and simplicity, show only vanilla DEXSeq.
ms$prefilter <- as.character(ms$prefilter)
pdf('analysis_sim_v3_perfHS.pdf')
print( ggplot(ms[species=='Human', ]) + gth3 +
facet_grid(type ~ tool, switch='y') +
geom_jitter(aes(x=effect, y=value, colour=reproducibility, shape=prefilter), width=0.1) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
coord_cartesian(ylim=c(0, 1)) +
scale_shape_manual(values=c(21, 4)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'black',
'Qrep >= 0.95 Rrep >= 0.85' = 'red',
'Qrep >= 0.8 Rrep >= 0.55' = 'blue')) +
labs(title='simulated human data', y='') +
guides(shape=guide_legend(nrow=2), colour=guide_legend(nrow=3)) +
theme(legend.position = "bottom") )
dev.off()## quartz_off_screen
## 2pdf('analysis_sim_v3_perfDS.pdf')
print( ggplot(ms[species=='Fruitfly', ]) + gth3 +
facet_grid(type ~ tool, switch='y') +
geom_jitter(aes(x=effect, y=value, colour=reproducibility, shape=prefilter), alpha=0.7, width=0.1) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
coord_cartesian(ylim=c(0, 1)) +
scale_shape_manual(values=c(21, 4)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'black',
'Qrep >= 0.95 Rrep >= 0.85' = 'red',
'Qrep >= 0.8 Rrep >= 0.55' = 'blue')) +
labs(title='simulated fruitfly data', y='') +
guides(shape=guide_legend(nrow=2), colour=guide_legend(nrow=3)) +
theme(legend.position = "bottom") )
dev.off()## quartz_off_screen
## 2# Human set quantified with Salmon.
hs <- ggplot(mega[quantification=='Salmon' & species=='Human', ]) + gth3 +
facet_grid(type ~ tool, switch='y') +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name="")) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated human dataset quantified with Salmon)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perfHS.gif", hs)
# Fly set quantified with Salmon.
ds <- ggplot(mega[quantification=='Salmon' & species=='Fruitfly', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name="")) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated fruitfly dataset quantified with Salmon)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perfDS.gif", ds)
# Human set quantified with Kallisto.
hk <- ggplot(mega[quantification=='Kallisto' & species=='Human', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name="")) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated human dataset quantified with Kallisto)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perfHK.gif", hk)
# Fruitfly set quantified with Kallisto.
dk <- ggplot(mega[quantification=='Kallisto' & species=='Fruitfly', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name="")) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated fruitfly dataset quantified with Kallisto)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perfDK.gif", dk)
- The sensitivity of RATs is not be affected by pre-filtering for these datasets, and the MCC is steady too. SUPPA2 and DEXSeq respond with a small reduction in sensitivity. In contrast DRIMSeq’s sensitivity improves noticeably with higher pre-filter values.
- SUPPA2 and RATs benefit the most by an FDR reduction from applying pre-filters.
- Post-filtering of effect size affects RATs and SUPPA2 in a similar way across all three metrics. DRIMSeq benefits from a small FDR reduction with the post-filter, but MCC and sensitivity change little.
- The stageR correction to DEXseq leads to a marginal improvement of MCC, with small reduction in FDR and sensitivity.
Rank
It is not easy to compare the points in the plots above, as many are quite close. In the tables below, we rank the tools, for each combination of prefilter and postfilter values.
For simpliciy, for RATs we are using the values from the strictest (ie. default) reproducibility thresholds when ranking (darkest red points in the plots). We are also limiting the ranking to the Salmon quantification, as Salmon and Kallisto resulted in practically indistinguishable preformance results.
As the ranks change easily with changes in the parameters, we use the sum of the ranks across the parameter combinations as an indication of which tool ranks high most consistently.
Human
Without prefilter
tallyh = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
sensh = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
fdrh = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
mcch = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Human') %>%
filter(prefilter == 0)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.612 D >= 0.2 Sensitivity
## 2 DEXSeq 0.570 D >= 0.2 Sensitivity
## 3 RATs(g) 0.563 D >= 0.2 Sensitivity
## 4 RATs(t) 0.549 D >= 0.2 Sensitivity
## 5 DRIMSeq(t) 0.513 D >= 0.2 Sensitivity
## 6 DRIMSeq(g) 0.505 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.021 D >= 0.2 FDR
## 2 DRIMSeq(t) 0.027 D >= 0.2 FDR
## 3 DEXSeq 0.032 D >= 0.2 FDR
## 4 RATs(g) 0.039 D >= 0.2 FDR
## 5 DRIMSeq(g) 0.062 D >= 0.2 FDR
## 6 SUPPA2 0.329 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.734 D >= 0.2 MCC
## 2 RATs(g) 0.726 D >= 0.2 MCC
## 3 RATs(t) 0.724 D >= 0.2 MCC
## 4 DRIMSeq(t) 0.691 D >= 0.2 MCC
## 5 DRIMSeq(g) 0.672 D >= 0.2 MCC
## 6 SUPPA2 0.623 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.788 D >= 0.1 Sensitivity
## 2 RATs(g) 0.716 D >= 0.1 Sensitivity
## 3 RATs(t) 0.692 D >= 0.1 Sensitivity
## 4 DEXSeq 0.599 D >= 0.1 Sensitivity
## 5 DRIMSeq(t) 0.558 D >= 0.1 Sensitivity
## 6 DRIMSeq(g) 0.555 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.061 D >= 0.1 FDR
## 2 DEXSeq 0.070 D >= 0.1 FDR
## 3 RATs(t) 0.071 D >= 0.1 FDR
## 4 DRIMSeq(g) 0.114 D >= 0.1 FDR
## 5 RATs(g) 0.282 D >= 0.1 FDR
## 6 SUPPA2 0.509 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.793 D >= 0.1 MCC
## 2 DEXSeq 0.737 D >= 0.1 MCC
## 3 DRIMSeq(t) 0.708 D >= 0.1 MCC
## 4 RATs(g) 0.703 D >= 0.1 MCC
## 5 DRIMSeq(g) 0.683 D >= 0.1 MCC
## 6 SUPPA2 0.598 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.871 D >= 0.05 Sensitivity
## 2 RATs(g) 0.864 D >= 0.05 Sensitivity
## 3 RATs(t) 0.777 D >= 0.05 Sensitivity
## 4 DEXSeq 0.602 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.574 D >= 0.05 Sensitivity
## 6 DRIMSeq(g) 0.566 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.103 D >= 0.05 FDR
## 2 DEXSeq 0.119 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.170 D >= 0.05 FDR
## 4 RATs(t) 0.306 D >= 0.05 FDR
## 5 SUPPA2 0.645 D >= 0.05 FDR
## 6 RATs(g) 0.670 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.720 D >= 0.05 MCC
## 2 DEXSeq 0.718 D >= 0.05 MCC
## 3 DRIMSeq(t) 0.699 D >= 0.05 MCC
## 4 DRIMSeq(g) 0.664 D >= 0.05 MCC
## 5 SUPPA2 0.524 D >= 0.05 MCC
## 6 RATs(g) 0.499 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}Prefilter > 10
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Human') %>%
filter(prefilter == 10)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.608 D >= 0.2 Sensitivity
## 2 DRIMSeq(g) 0.596 D >= 0.2 Sensitivity
## 3 DRIMSeq(t) 0.589 D >= 0.2 Sensitivity
## 4 DEXSeq 0.567 D >= 0.2 Sensitivity
## 5 RATs(g) 0.563 D >= 0.2 Sensitivity
## 6 RATs(t) 0.549 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.014 D >= 0.2 FDR
## 2 RATs(g) 0.016 D >= 0.2 FDR
## 3 DEXSeq 0.027 D >= 0.2 FDR
## 4 SUPPA2 0.046 D >= 0.2 FDR
## 5 DRIMSeq(t) 0.050 D >= 0.2 FDR
## 6 DRIMSeq(g) 0.071 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.750 D >= 0.2 MCC
## 2 RATs(g) 0.736 D >= 0.2 MCC
## 3 DEXSeq 0.734 D >= 0.2 MCC
## 4 DRIMSeq(t) 0.729 D >= 0.2 MCC
## 5 RATs(t) 0.727 D >= 0.2 MCC
## 6 DRIMSeq(g) 0.725 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.787 D >= 0.1 Sensitivity
## 2 RATs(g) 0.716 D >= 0.1 Sensitivity
## 3 RATs(t) 0.692 D >= 0.1 Sensitivity
## 4 DRIMSeq(g) 0.651 D >= 0.1 Sensitivity
## 5 DRIMSeq(t) 0.641 D >= 0.1 Sensitivity
## 6 DEXSeq 0.599 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.058 D >= 0.1 FDR
## 2 RATs(t) 0.061 D >= 0.1 FDR
## 3 DRIMSeq(t) 0.109 D >= 0.1 FDR
## 4 DRIMSeq(g) 0.177 D >= 0.1 FDR
## 5 RATs(g) 0.211 D >= 0.1 FDR
## 6 SUPPA2 0.255 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.798 D >= 0.1 MCC
## 2 SUPPA2 0.750 D >= 0.1 MCC
## 3 DEXSeq 0.742 D >= 0.1 MCC
## 4 RATs(g) 0.739 D >= 0.1 MCC
## 5 DRIMSeq(t) 0.735 D >= 0.1 MCC
## 6 DRIMSeq(g) 0.707 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.878 D >= 0.05 Sensitivity
## 2 RATs(g) 0.864 D >= 0.05 Sensitivity
## 3 RATs(t) 0.778 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.656 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.649 D >= 0.05 Sensitivity
## 6 DEXSeq 0.601 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.106 D >= 0.05 FDR
## 2 DRIMSeq(t) 0.170 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.266 D >= 0.05 FDR
## 4 RATs(t) 0.303 D >= 0.05 FDR
## 5 SUPPA2 0.561 D >= 0.05 FDR
## 6 RATs(g) 0.655 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.723 D >= 0.05 MCC
## 2 RATs(t) 0.722 D >= 0.05 MCC
## 3 DRIMSeq(t) 0.710 D >= 0.05 MCC
## 4 DRIMSeq(g) 0.663 D >= 0.05 MCC
## 5 SUPPA2 0.589 D >= 0.05 MCC
## 6 RATs(g) 0.513 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}Prefilter > 50
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Human') %>%
filter(prefilter == 50)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.632 D >= 0.2 Sensitivity
## 2 DRIMSeq(t) 0.616 D >= 0.2 Sensitivity
## 3 SUPPA2 0.600 D >= 0.2 Sensitivity
## 4 RATs(g) 0.562 D >= 0.2 Sensitivity
## 5 RATs(t) 0.549 D >= 0.2 Sensitivity
## 6 DEXSeq 0.523 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.004 D >= 0.2 FDR
## 2 RATs(t) 0.004 D >= 0.2 FDR
## 3 DRIMSeq(t) 0.033 D >= 0.2 FDR
## 4 DRIMSeq(g) 0.046 D >= 0.2 FDR
## 5 SUPPA2 0.053 D >= 0.2 FDR
## 6 DEXSeq 0.053 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.754 D >= 0.2 MCC
## 2 DRIMSeq(t) 0.749 D >= 0.2 MCC
## 3 SUPPA2 0.748 D >= 0.2 MCC
## 4 RATs(g) 0.740 D >= 0.2 MCC
## 5 RATs(t) 0.731 D >= 0.2 MCC
## 6 DEXSeq 0.694 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.776 D >= 0.1 Sensitivity
## 2 RATs(g) 0.715 D >= 0.1 Sensitivity
## 3 RATs(t) 0.692 D >= 0.1 Sensitivity
## 4 DRIMSeq(g) 0.678 D >= 0.1 Sensitivity
## 5 DRIMSeq(t) 0.665 D >= 0.1 Sensitivity
## 6 DEXSeq 0.552 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.036 D >= 0.1 FDR
## 2 RATs(g) 0.073 D >= 0.1 FDR
## 3 DRIMSeq(t) 0.081 D >= 0.1 FDR
## 4 DEXSeq 0.089 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.126 D >= 0.1 FDR
## 6 SUPPA2 0.149 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.809 D >= 0.1 MCC
## 2 SUPPA2 0.807 D >= 0.1 MCC
## 3 RATs(g) 0.806 D >= 0.1 MCC
## 4 DRIMSeq(t) 0.758 D >= 0.1 MCC
## 5 DRIMSeq(g) 0.743 D >= 0.1 MCC
## 6 DEXSeq 0.698 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.902 D >= 0.05 Sensitivity
## 2 RATs(g) 0.864 D >= 0.05 Sensitivity
## 3 RATs(t) 0.777 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.683 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.671 D >= 0.05 Sensitivity
## 6 DEXSeq 0.555 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.122 D >= 0.05 FDR
## 2 DEXSeq 0.146 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.178 D >= 0.05 FDR
## 4 RATs(t) 0.306 D >= 0.05 FDR
## 5 SUPPA2 0.562 D >= 0.05 FDR
## 6 RATs(g) 0.670 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.740 D >= 0.05 MCC
## 2 RATs(t) 0.720 D >= 0.05 MCC
## 3 DRIMSeq(g) 0.718 D >= 0.05 MCC
## 4 DEXSeq 0.676 D >= 0.05 MCC
## 5 SUPPA2 0.613 D >= 0.05 MCC
## 6 RATs(g) 0.499 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}Cumulative rank
print( data.frame(Overall = tallyh,
Sensitivity = sensh,
FDR = fdrh,
MCC = mcch) )## Overall Sensitivity FDR MCC
## RATs(t) 76 33 22 21
## RATs(g) 98 24 37 37
## DRIMSeq(t) 90 40 21 29
## DRIMSeq(g) 113 37 37 39
## SUPPA2 94 11 48 35
## DEXSeq 96 44 24 28Although the tool occupying the top rank changes easily with the filtering parameters of the comparison, a crude sum of the ranks from the previous tables across the comparisons shows that overall the transcript-level method in RATs consistently ranks higher than the other tools, even when it is not in the top spot.
Fruitfly
Without prefilter
tallyd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
sensd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
fdrd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
mccd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Fruitfly') %>%
filter(prefilter == 0)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.615 D >= 0.2 Sensitivity
## 2 DRIMSeq(g) 0.601 D >= 0.2 Sensitivity
## 3 DRIMSeq(t) 0.595 D >= 0.2 Sensitivity
## 4 DEXSeq 0.577 D >= 0.2 Sensitivity
## 5 RATs(g) 0.535 D >= 0.2 Sensitivity
## 6 RATs(t) 0.526 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.009 D >= 0.2 FDR
## 2 RATs(g) 0.015 D >= 0.2 FDR
## 3 DEXSeq 0.015 D >= 0.2 FDR
## 4 DRIMSeq(t) 0.029 D >= 0.2 FDR
## 5 DRIMSeq(g) 0.036 D >= 0.2 FDR
## 6 SUPPA2 0.128 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.741 D >= 0.2 MCC
## 2 DRIMSeq(g) 0.721 D >= 0.2 MCC
## 3 DRIMSeq(t) 0.721 D >= 0.2 MCC
## 4 SUPPA2 0.716 D >= 0.2 MCC
## 5 RATs(g) 0.713 D >= 0.2 MCC
## 6 RATs(t) 0.709 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.765 D >= 0.1 Sensitivity
## 2 RATs(g) 0.666 D >= 0.1 Sensitivity
## 3 RATs(t) 0.653 D >= 0.1 Sensitivity
## 4 DRIMSeq(t) 0.649 D >= 0.1 Sensitivity
## 5 DRIMSeq(g) 0.640 D >= 0.1 Sensitivity
## 6 DEXSeq 0.616 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.024 D >= 0.1 FDR
## 2 DEXSeq 0.033 D >= 0.1 FDR
## 3 RATs(g) 0.043 D >= 0.1 FDR
## 4 DRIMSeq(t) 0.056 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.071 D >= 0.1 FDR
## 6 SUPPA2 0.223 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.787 D >= 0.1 MCC
## 2 RATs(g) 0.786 D >= 0.1 MCC
## 3 DEXSeq 0.759 D >= 0.1 MCC
## 4 SUPPA2 0.753 D >= 0.1 MCC
## 5 DRIMSeq(t) 0.743 D >= 0.1 MCC
## 6 DRIMSeq(g) 0.729 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.801 D >= 0.05 Sensitivity
## 2 SUPPA2 0.786 D >= 0.05 Sensitivity
## 3 RATs(t) 0.754 D >= 0.05 Sensitivity
## 4 DRIMSeq(t) 0.652 D >= 0.05 Sensitivity
## 5 DRIMSeq(g) 0.642 D >= 0.05 Sensitivity
## 6 DEXSeq 0.621 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.052 D >= 0.05 FDR
## 2 DRIMSeq(t) 0.072 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.099 D >= 0.05 FDR
## 4 RATs(t) 0.118 D >= 0.05 FDR
## 5 RATs(g) 0.222 D >= 0.05 FDR
## 6 SUPPA2 0.248 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.803 D >= 0.05 MCC
## 2 RATs(g) 0.773 D >= 0.05 MCC
## 3 DEXSeq 0.754 D >= 0.05 MCC
## 4 SUPPA2 0.750 D >= 0.05 MCC
## 5 DRIMSeq(t) 0.736 D >= 0.05 MCC
## 6 DRIMSeq(g) 0.715 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}Prefilter > 10
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Fruitfly') %>%
filter(prefilter == 10)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.635 D >= 0.2 Sensitivity
## 2 SUPPA2 0.618 D >= 0.2 Sensitivity
## 3 DRIMSeq(t) 0.614 D >= 0.2 Sensitivity
## 4 DEXSeq 0.588 D >= 0.2 Sensitivity
## 5 RATs(g) 0.535 D >= 0.2 Sensitivity
## 6 RATs(t) 0.526 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.008 D >= 0.2 FDR
## 2 RATs(g) 0.009 D >= 0.2 FDR
## 3 DEXSeq 0.018 D >= 0.2 FDR
## 4 SUPPA2 0.026 D >= 0.2 FDR
## 5 DRIMSeq(t) 0.029 D >= 0.2 FDR
## 6 DRIMSeq(g) 0.033 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.759 D >= 0.2 MCC
## 2 DEXSeq 0.747 D >= 0.2 MCC
## 3 DRIMSeq(g) 0.740 D >= 0.2 MCC
## 4 DRIMSeq(t) 0.727 D >= 0.2 MCC
## 5 RATs(g) 0.715 D >= 0.2 MCC
## 6 RATs(t) 0.709 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.772 D >= 0.1 Sensitivity
## 2 DRIMSeq(g) 0.683 D >= 0.1 Sensitivity
## 3 DRIMSeq(t) 0.682 D >= 0.1 Sensitivity
## 4 RATs(g) 0.666 D >= 0.1 Sensitivity
## 5 DEXSeq 0.663 D >= 0.1 Sensitivity
## 6 RATs(t) 0.653 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.021 D >= 0.1 FDR
## 2 RATs(g) 0.036 D >= 0.1 FDR
## 3 DEXSeq 0.058 D >= 0.1 FDR
## 4 DRIMSeq(t) 0.066 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.073 D >= 0.1 FDR
## 6 SUPPA2 0.136 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.798 D >= 0.1 MCC
## 2 RATs(g) 0.789 D >= 0.1 MCC
## 3 RATs(t) 0.788 D >= 0.1 MCC
## 4 DEXSeq 0.778 D >= 0.1 MCC
## 5 DRIMSeq(t) 0.753 D >= 0.1 MCC
## 6 DRIMSeq(g) 0.749 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.801 D >= 0.05 Sensitivity
## 2 SUPPA2 0.795 D >= 0.05 Sensitivity
## 3 RATs(t) 0.754 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.685 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.684 D >= 0.05 Sensitivity
## 6 DEXSeq 0.670 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.086 D >= 0.05 FDR
## 2 DEXSeq 0.088 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.098 D >= 0.05 FDR
## 4 RATs(t) 0.116 D >= 0.05 FDR
## 5 SUPPA2 0.171 D >= 0.05 FDR
## 6 RATs(g) 0.218 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.804 D >= 0.05 MCC
## 2 SUPPA2 0.792 D >= 0.05 MCC
## 3 RATs(g) 0.775 D >= 0.05 MCC
## 4 DEXSeq 0.768 D >= 0.05 MCC
## 5 DRIMSeq(t) 0.742 D >= 0.05 MCC
## 6 DRIMSeq(g) 0.736 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}Prefilter > 50
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Fruitfly') %>%
filter(prefilter == 50)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.647 D >= 0.2 Sensitivity
## 2 DRIMSeq(t) 0.631 D >= 0.2 Sensitivity
## 3 SUPPA2 0.628 D >= 0.2 Sensitivity
## 4 DEXSeq 0.556 D >= 0.2 Sensitivity
## 5 RATs(g) 0.535 D >= 0.2 Sensitivity
## 6 RATs(t) 0.526 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.007 D >= 0.2 FDR
## 2 RATs(t) 0.008 D >= 0.2 FDR
## 3 SUPPA2 0.012 D >= 0.2 FDR
## 4 DRIMSeq(t) 0.024 D >= 0.2 FDR
## 5 DRIMSeq(g) 0.026 D >= 0.2 FDR
## 6 DEXSeq 0.030 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.773 D >= 0.2 MCC
## 2 DRIMSeq(g) 0.738 D >= 0.2 MCC
## 3 DRIMSeq(t) 0.728 D >= 0.2 MCC
## 4 DEXSeq 0.721 D >= 0.2 MCC
## 5 RATs(g) 0.716 D >= 0.2 MCC
## 6 RATs(t) 0.709 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.781 D >= 0.1 Sensitivity
## 2 DRIMSeq(g) 0.703 D >= 0.1 Sensitivity
## 3 DRIMSeq(t) 0.700 D >= 0.1 Sensitivity
## 4 RATs(g) 0.666 D >= 0.1 Sensitivity
## 5 RATs(t) 0.653 D >= 0.1 Sensitivity
## 6 DEXSeq 0.612 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.017 D >= 0.1 FDR
## 2 RATs(g) 0.023 D >= 0.1 FDR
## 3 DRIMSeq(t) 0.054 D >= 0.1 FDR
## 4 DEXSeq 0.058 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.061 D >= 0.1 FDR
## 6 SUPPA2 0.070 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.839 D >= 0.1 MCC
## 2 RATs(g) 0.795 D >= 0.1 MCC
## 3 RATs(t) 0.790 D >= 0.1 MCC
## 4 DRIMSeq(t) 0.757 D >= 0.1 MCC
## 5 DRIMSeq(g) 0.756 D >= 0.1 MCC
## 6 DEXSeq 0.745 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.817 D >= 0.05 Sensitivity
## 2 RATs(g) 0.801 D >= 0.05 Sensitivity
## 3 RATs(t) 0.754 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.704 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.702 D >= 0.05 Sensitivity
## 6 DEXSeq 0.616 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.063 D >= 0.05 FDR
## 2 DRIMSeq(g) 0.076 D >= 0.05 FDR
## 3 DEXSeq 0.083 D >= 0.05 FDR
## 4 RATs(t) 0.118 D >= 0.05 FDR
## 5 SUPPA2 0.149 D >= 0.05 FDR
## 6 RATs(g) 0.222 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.818 D >= 0.05 MCC
## 2 RATs(t) 0.803 D >= 0.05 MCC
## 3 RATs(g) 0.773 D >= 0.05 MCC
## 4 DRIMSeq(t) 0.753 D >= 0.05 MCC
## 5 DRIMSeq(g) 0.746 D >= 0.05 MCC
## 6 DEXSeq 0.737 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}Cumulative rank
print( data.frame(Overall = tallyd,
Sensitivity = sensd,
FDR = fdrd,
MCC = mccd) )## Overall Sensitivity FDR MCC
## RATs(t) 89 41 19 29
## RATs(g) 87 29 29 29
## DRIMSeq(t) 98 32 28 38
## DRIMSeq(g) 106 26 39 41
## SUPPA2 80 14 47 19
## DEXSeq 107 47 27 33Both datasets
tmp <- data.frame(Cumulative_Rank= c(tallyh + tallyd, sensh + sensd, fdrh + fdrd, mcch + mccd),
Type=factor(rep(c('Overall', 'Sensitivity', 'FDR', 'MCC'), each=6), levels=c('Overall', 'Sensitivity', 'FDR', 'MCC')),
Tool=c(names(tallyh), names(sensh), names(fdrh), names(mcch)) )
p <-ggplot(tmp, aes(x=Tool, fill=Type, y=Cumulative_Rank)) + gth +
geom_bar(stat='identity', position='dodge') +
scale_y_reverse() +
scale_x_discrete(position = "top") +
scale_fill_manual(values=c(Overall='grey70', Sensitivity='orange', FDR='red', MCC='blue'))
print(p)
pdf('cumrank.pdf')
print(p)
dev.off()## quartz_off_screen
## 2Comparison between the tools at given thresholds, accounting for genes excluded from the simulations.
The simulated sets were focused on coding genes only. Non-coding genes would be expected to appear as non-expressed and are added to the FP and TN.
Assistant functions
Redefine the calculations of TP,FP,etc to include non-coding genes as expected Negatives.
# Redefine the metrics function
measureup <- function(x, tool, fx, qt=0.95, rt=0.85, pre=0) {
df <- countup(x, tool, fx, qt, rt)
dt <- data.table(value= c(round(df[, tp/(tp+fn)], digits=3),
round(df[, (fp+fna)/(tp+fp+fna)], digits=3),
round(df[, as.double(tp*(tn+tna) - (fp+fna)*fn) / sqrt( as.double(tp + fp+fna)*as.double(tp + fn)*as.double(tn+tna + fp+fna)*as.double(tn+tna + fn) )], digits=3)),
type= factor(c(rep('Sensitivity', 4), rep('FDR', 4), rep('MCC', 4)), levels=c('Sensitivity', 'FDR', 'MCC')),
tool= factor(rep(tool, 12), levels=c('DRIMSeq(g)', 'DRIMSeq(t)', 'RATs(g)', 'RATs(t)', 'SUPPA2', 'DEXSeq')),
species=rep(df$set, 3),
quantification=rep(df$quant, 3),
effect=rep(paste('D >= ', fx), 12) ,
reproducibility=factor(rep(paste('Qrep', ifelse(is.na(qt), '', '>='), qt, 'Rrep', ifelse(is.na(rt), '', '>='), rt), 12)),
prefilter=pre)
if (tool=='DEXSeq')
dt[, tool := rep(c('DEXSeq', 'DEXsgr'), 6)]
return(dt)
}Performance at pval < 0.05
RATs
# RATs Dprop >= 0.20
rgd <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.95, rt=0.85, pre=0)
rtd <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.95, rt=0.85, pre=0)
# more permissive reproducibility
rgd1 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.90, rt=0.75, pre=0)
rtd1 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.90, rt=0.75, pre=0)
rgd2 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.85, rt=0.65, pre=0)
rtd2 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.85, rt=0.65, pre=0)
rgd3 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(g)', fx=0.20, qt=0.80, rt=0.55, pre=0)
rtd3 <- measureup(list(rdk,rds,rhk,rhs), 'RATs(t)', fx=0.20, qt=0.80, rt=0.55, pre=0)
# RATs Dprop >= 0.10
rgf <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.95, rt=0.85, pre=0)
rtf <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.95, rt=0.85, pre=0)
# more permissive reproducibility
rgf1 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.90, rt=0.75, pre=0)
rtf1 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.90, rt=0.75, pre=0)
rgf2 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.85, rt=0.65, pre=0)
rtf2 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.85, rt=0.65, pre=0)
rgf3 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(g)', fx=0.10, qt=0.80, rt=0.55, pre=0)
rtf3 <- measureup(list(rdkth,rdsth,rhkth,rhsth), 'RATs(t)', fx=0.10, qt=0.80, rt=0.55, pre=0)
# RATs Dprop >= 0.05
rgff <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.95, rt=0.85, pre=0)
rtff <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.95, rt=0.85, pre=0)
# more permissive reproducibility
rgff1 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.90, rt=0.75, pre=0)
rtff1 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.90, rt=0.75, pre=0)
rgff2 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.85, rt=0.65, pre=0)
rtff2 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.85, rt=0.65, pre=0)
rgff3 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.80, rt=0.55, pre=0)
rtff3 <- measureup(list(rdkth2,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.80, rt=0.55, pre=0)# RATs Dprop >= 0.20
rgd10 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.95, rt=0.85, pre=10)
rtd10 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.95, rt=0.85, pre=10)
# more permissive reproducibility
rgd110 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.90, rt=0.75, pre=10)
rtd110 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.90, rt=0.75, pre=10)
rgd210 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.85, rt=0.65, pre=10)
rtd210 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.85, rt=0.65, pre=10)
rgd310 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(g)', fx=0.20, qt=0.80, rt=0.55, pre=10)
rtd310 <- measureup(list(rdk10,rds10,rhk10,rhs10), 'RATs(t)', fx=0.20, qt=0.80, rt=0.55, pre=10)
# RATs Dprop >= 0.10
rgf10 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.95, rt=0.85, pre=10)
rtf10 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.95, rt=0.85, pre=10)
# more permissive reproducibility
rgf110 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.90, rt=0.75, pre=10)
rtf110 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.90, rt=0.75, pre=10)
rgf210 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.85, rt=0.65, pre=10)
rtf210 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.85, rt=0.65, pre=10)
rgf310 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(g)', fx=0.10, qt=0.80, rt=0.55, pre=10)
rtf310 <- measureup(list(rdkth10,rdsth10,rhkth10,rhsth10), 'RATs(t)', fx=0.10, qt=0.80, rt=0.55, pre=10)
# RATs Dprop >= 0.05
rgff10 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.95, rt=0.85, pre=10)
rtff10 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.95, rt=0.85, pre=10)
# more permissive reproducibility
rgff110 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.90, rt=0.75, pre=10)
rtff110 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.90, rt=0.75, pre=10)
rgff210 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.85, rt=0.65, pre=10)
rtff210 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.85, rt=0.65, pre=10)
rgff310 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(g)', fx=0.05, qt=0.80, rt=0.55, pre=10)
rtff310 <- measureup(list(rdkth210,rdsth210,rhkth210,rhsth210), 'RATs(t)', fx=0.05, qt=0.80, rt=0.55, pre=10)# RATs Dprop >= 0.20
rgd50 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.95, rt=0.85, pre=50)
rtd50 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.95, rt=0.85, pre=50)
# more permissive reproducibility
rgd150 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.90, rt=0.75, pre=50)
rtd150 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.90, rt=0.75, pre=50)
rgd250 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.85, rt=0.65, pre=50)
rtd250 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.85, rt=0.65, pre=50)
rgd350 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(g)', fx=0.20, qt=0.80, rt=0.55, pre=50)
rtd350 <- measureup(list(rdk50,rds50,rhk50,rhs50), 'RATs(t)', fx=0.20, qt=0.80, rt=0.55, pre=50)
# RATs Dprop >= 0.10
rgf50 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.95, rt=0.85, pre=50)
rtf50 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.95, rt=0.85, pre=50)
# more permissive reproducibility
rgf150 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.90, rt=0.75, pre=50)
rtf150 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.90, rt=0.75, pre=50)
rgf250 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.85, rt=0.65, pre=50)
rtf250 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.85, rt=0.65, pre=50)
rgf350 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(g)', fx=0.10, qt=0.80, rt=0.55, pre=50)
rtf350 <- measureup(list(rdkth50,rdsth50,rhkth50,rhsth50), 'RATs(t)', fx=0.10, qt=0.80, rt=0.55, pre=50)
# RATs Dprop >= 0.05
rgff50 <- measureup(list(rdkth250,rdsth2,rhkth2,rhsth2), 'RATs(g)', fx=0.05, qt=0.95, rt=0.85, pre=50)
rtff50 <- measureup(list(rdkth250,rdsth2,rhkth2,rhsth2), 'RATs(t)', fx=0.05, qt=0.95, rt=0.85, pre=50)
# more permissive reproducibility
rgff150 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(g)', fx=0.05, qt=0.90, rt=0.75, pre=50)
rtff150 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(t)', fx=0.05, qt=0.90, rt=0.75, pre=50)
rgff250 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(g)', fx=0.05, qt=0.85, rt=0.65, pre=50)
rtff250 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(t)', fx=0.05, qt=0.85, rt=0.65, pre=50)
rgff350 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(g)', fx=0.05, qt=0.80, rt=0.55, pre=50)
rtff350 <- measureup(list(rdkth250,rdsth250,rhkth250,rhsth250), 'RATs(t)', fx=0.05, qt=0.80, rt=0.55, pre=50)SUPPA2
# SUPPA2 does not bootstrap the reproducibility.
# SUPPA2 dPSI >= 0.20
sd <- measureup(list(sdk,sds,shk,shs), 'SUPPA2', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0)
# SUPPA2 dPSI >= 0.10
sf <- measureup(list(sdk,sds,shk,shs), 'SUPPA2', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0)
# SUPPA2 dPSI >= 0.05
sff <- measureup(list(sdk,sds,shk,shs), 'SUPPA2', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0)# SUPPA2 does not bootstrap the reproducibility.
# SUPPA2 dPSI >= 0.20
sd10 <- measureup(list(sdk10,sds10,shk10,shs10), 'SUPPA2', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10)
# SUPPA2 dPSI >= 0.10
sf10 <- measureup(list(sdk10,sds10,shk10,shs10), 'SUPPA2', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10)
# SUPPA2 dPSI >= 0.05
sff10 <- measureup(list(sdk10,sds10,shk10,shs10), 'SUPPA2', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10)# SUPPA2 does not bootstrap the reproducibility.
# SUPPA2 dPSI >= 0.20
sd50 <- measureup(list(sdk50,sds50,shk50,shs50), 'SUPPA2', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50)
# SUPPA2 dPSI >= 0.10
sf50 <- measureup(list(sdk50,sds50,shk50,shs50), 'SUPPA2', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50)
# SUPPA2 dPSI >= 0.05
sff50 <- measureup(list(sdk50,sds50,shk50,shs50), 'SUPPA2', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50)DRIMSeq
# DRIMSeq does not bootstrap the reproducibility.
# Gene-level
# Dprop >= 0.20
dd <- measureup(list(ddk,dds,dhk,dhs), 'DRIMSeq(g)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.10
df <- measureup(list(ddk,dds,dhk,dhs), 'DRIMSeq(g)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.05
dff <- measureup(list(ddk,dds,dhk,dhs), 'DRIMSeq(g)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0)
# Transcript-level
# Dprop >= 0.20
ddt <- measureup(list(ddkt,ddst,dhkt,dhst), 'DRIMSeq(t)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.10
dft <- measureup(list(ddkt,ddst,dhkt,dhst), 'DRIMSeq(t)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0)
# Dprop >= 0.05
dfft <- measureup(list(ddkt,ddst,dhkt,dhst), 'DRIMSeq(t)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0)# DRIMSeq does not bootstrap the reproducibility.
# Gene-level
# Dprop >= 0.20
dd10 <- measureup(list(ddk10,dds10,dhk10,dhs10), 'DRIMSeq(g)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.10
df10 <- measureup(list(ddk10,dds10,dhk10,dhs10), 'DRIMSeq(g)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.05
dff10 <- measureup(list(ddk10,dds10,dhk10,dhs10), 'DRIMSeq(g)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10)
# Transcript-level
# Dprop >= 0.20
ddt10 <- measureup(list(ddkt10,ddst10,dhkt10,dhst10), 'DRIMSeq(t)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.10
dft10 <- measureup(list(ddkt10,ddst10,dhkt10,dhst10), 'DRIMSeq(t)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10)
# Dprop >= 0.05
dfft10 <- measureup(list(ddkt10,ddst10,dhkt10,dhst10), 'DRIMSeq(t)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10)# DRIMSeq does not bootstrap the reproducibility.
# Gene-level
# Dprop >= 0.20
dd50 <- measureup(list(ddk50,dds50,dhk50,dhs50), 'DRIMSeq(g)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.10
df50 <- measureup(list(ddk50,dds50,dhk50,dhs50), 'DRIMSeq(g)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.05
dff50 <- measureup(list(ddk50,dds50,dhk50,dhs50), 'DRIMSeq(g)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50)
# Transcript-level
# Dprop >= 0.20
ddt50 <- measureup(list(ddkt50,ddst50,dhkt50,dhst50), 'DRIMSeq(t)', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.10
dft50 <- measureup(list(ddkt50,ddst50,dhkt50,dhst50), 'DRIMSeq(t)', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50)
# Dprop >= 0.05
dfft50 <- measureup(list(ddkt50,ddst50,dhkt50,dhst50), 'DRIMSeq(t)', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50)DEXSeq
# DEXSeq does not bootstrap the reproducibility.
# DEXSeq with and without stageR corrections. Dprop >= 0.20
xd <- measureup(list(xds, xsds, xhs, xshs), 'DEXSeq', fx=0.20, qt=NA_real_, rt=NA_real_, pre=0 )
# DEXSeq with and without stageR corrections. Dprop >= 0.10
xf <- measureup(list(xds, xsds, xhs, xshs), 'DEXSeq', fx=0.10, qt=NA_real_, rt=NA_real_, pre=0 )
# DEXSeq with and without stageR corrections. Dprop >= 0.05
xff <- measureup(list(xds, xsds, xhs, xshs), 'DEXSeq', fx=0.05, qt=NA_real_, rt=NA_real_, pre=0 )# DEXSeq does not bootstrap the reproducibility.
# DEXSeq with and without stageR corrections. Dprop >= 0.20
xd10 <- measureup(list(xds10, xsds10, xhs10, xshs10), 'DEXSeq', fx=0.20, qt=NA_real_, rt=NA_real_, pre=10 )
# DEXSeq with and without stageR corrections. Dprop >= 0.10
xf10 <- measureup(list(xds10, xsds10, xhs10, xshs10), 'DEXSeq', fx=0.10, qt=NA_real_, rt=NA_real_, pre=10 )
# DEXSeq with and without stageR corrections. Dprop >= 0.05
xff10 <- measureup(list(xds10, xsds10, xhs10, xshs10), 'DEXSeq', fx=0.05, qt=NA_real_, rt=NA_real_, pre=10 )# DEXSeq does not bootstrap the reproducibility.
# DEXSeq with and without stageR corrections. Dprop >= 0.20
xd50 <- measureup(list(xds50, xsds50, xhs50, xshs50), 'DEXSeq', fx=0.20, qt=NA_real_, rt=NA_real_, pre=50 )
# DEXSeq with and without stageR corrections. Dprop >= 0.10
xf50 <- measureup(list(xds50, xsds50, xhs50, xshs50), 'DEXSeq', fx=0.10, qt=NA_real_, rt=NA_real_, pre=50 )
# DEXSeq with and without stageR corrections. Dprop >= 0.05
xff50 <- measureup(list(xds50, xsds50, xhs50, xshs50), 'DEXSeq', fx=0.05, qt=NA_real_, rt=NA_real_, pre=50 )Plot
# Easier to work with ggplot2 if everything is in one table.
mega <- rbindlist(list(rgd, rtd, rgd1, rtd1, rgd2, rtd2, rgd3, rtd3,
rgf, rtf, rgf1, rtf1, rgf2, rtf2, rgf3, rtf3,
rgff, rtff, rgff1, rtff1, rgff2, rtff2, rgff3, rtff3,
sd, sf, sff, dd, df, dff, ddt, dft, dfft, xd, xf, xff,
rgd10, rtd10, rgd110, rtd110, rgd210, rtd210, rgd310, rtd310,
rgf10, rtf10, rgf110, rtf110, rgf210, rtf210, rgf310, rtf310,
rgff10, rtff10, rgff110, rtff110, rgff210, rtff210, rgff310, rtff310,
sd10, sf10, sff10, dd10, df10, dff10, ddt10, dft10, dfft10, xd10, xf10, xff10,
rgd50, rtd50, rgd150, rtd150, rgd250, rtd250, rgd350, rtd350,
rgf50, rtf50, rgf150, rtf150, rgf250, rtf250, rgf350, rtf350,
rgff50, rtff50, rgff150, rtff150, rgff250, rtff250, rgff350, rtff350,
sd50, sf50, sff50, dd50, df50, dff50, ddt50, dft50, dfft50, xd50, xf50, xff50))Plot.
# Human set quantified with Salmon.
hs <- ggplot(mega[quantification=='Salmon' & species=='Human', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name='')) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated human dataset quantified with Salmon)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perf+HS.gif", hs)
hs2 <- ggplot(mega[quantification=='Salmon' & species=='Human' & tool != 'DEXsgr', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name='')) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated human dataset quantified with Salmon)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perf+HS-.gif", hs2)
# Fly set quantified with Salmon.
ds <- ggplot(mega[quantification=='Salmon' & species=='Fruitfly', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name='')) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated fruitfly dataset quantified with Salmon)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perf+DS.gif", ds)
ds2 <- ggplot(mega[quantification=='Salmon' & species=='Fruitfly' & tool != 'DEXsgr', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name='')) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated fruitfly dataset quantified with Salmon)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perf+DS-.gif", ds2)
# Human set quantified with Kallisto.
hk <- ggplot(mega[quantification=='Kallisto' & species=='Human', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name='')) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated human dataset quantified with Kallisto)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perf+HK.gif", hk)
# Fruitfly set quantified with Kallisto.
dk <- ggplot(mega[quantification=='Kallisto' & species=='Fruitfly', ]) + gth3 +
facet_grid(type ~ tool, switch='y', scale="free_y") +
transition_states(prefilter, transition_length = 0.5, state_length = 3) +
geom_point(aes(x=effect, y=value, colour=reproducibility, shape=tool), shape=20) +
scale_y_continuous(expand=c(0,0), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1), sec.axis=dup_axis(name='')) +
coord_cartesian(ylim=c(0, 1)) +
scale_colour_manual(values=c('Qrep NA Rrep NA' = 'blue',
'Qrep >= 0.95 Rrep >= 0.85' = 'darkred',
'Qrep >= 0.9 Rrep >= 0.75' = 'red',
'Qrep >= 0.85 Rrep >= 0.65' = 'orange',
'Qrep >= 0.8 Rrep >= 0.55' = 'gold')) +
labs(title='DTU performance at typical significance level of 0.05', subtitle='(simulated fruitfly dataset quantified with Kallisto)', y="Abundance prefilter : {closest_state}") +
guides(colour=guide_legend(nrow=2)) +
theme(legend.position = "bottom")
anim_save("analysis_sim_v3_perf+DK.gif", dk)
The inclusion of the excluded genes as additional negative cases does not alter the result significantly.
Rank
Human
Without prefilter
tallyh = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
sensh = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
fdrh = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
mcch = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Human') %>%
filter(prefilter == 0)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.612 D >= 0.2 Sensitivity
## 2 DEXSeq 0.570 D >= 0.2 Sensitivity
## 3 RATs(g) 0.563 D >= 0.2 Sensitivity
## 4 RATs(t) 0.549 D >= 0.2 Sensitivity
## 5 DRIMSeq(t) 0.513 D >= 0.2 Sensitivity
## 6 DRIMSeq(g) 0.505 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.028 D >= 0.2 FDR
## 2 DRIMSeq(t) 0.033 D >= 0.2 FDR
## 3 DEXSeq 0.039 D >= 0.2 FDR
## 4 RATs(g) 0.049 D >= 0.2 FDR
## 5 DRIMSeq(g) 0.068 D >= 0.2 FDR
## 6 SUPPA2 0.358 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.737 D >= 0.2 MCC
## 2 RATs(g) 0.727 D >= 0.2 MCC
## 3 RATs(t) 0.727 D >= 0.2 MCC
## 4 DRIMSeq(t) 0.689 D >= 0.2 MCC
## 5 DRIMSeq(g) 0.678 D >= 0.2 MCC
## 6 SUPPA2 0.621 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.788 D >= 0.1 Sensitivity
## 2 RATs(g) 0.716 D >= 0.1 Sensitivity
## 3 RATs(t) 0.692 D >= 0.1 Sensitivity
## 4 DEXSeq 0.599 D >= 0.1 Sensitivity
## 5 DRIMSeq(t) 0.558 D >= 0.1 Sensitivity
## 6 DRIMSeq(g) 0.555 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.067 D >= 0.1 FDR
## 2 DEXSeq 0.076 D >= 0.1 FDR
## 3 RATs(t) 0.077 D >= 0.1 FDR
## 4 DRIMSeq(g) 0.120 D >= 0.1 FDR
## 5 RATs(g) 0.291 D >= 0.1 FDR
## 6 SUPPA2 0.525 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.796 D >= 0.1 MCC
## 2 DEXSeq 0.741 D >= 0.1 MCC
## 3 DRIMSeq(t) 0.706 D >= 0.1 MCC
## 4 RATs(g) 0.705 D >= 0.1 MCC
## 5 DRIMSeq(g) 0.690 D >= 0.1 MCC
## 6 SUPPA2 0.603 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.871 D >= 0.05 Sensitivity
## 2 RATs(g) 0.864 D >= 0.05 Sensitivity
## 3 RATs(t) 0.777 D >= 0.05 Sensitivity
## 4 DEXSeq 0.602 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.574 D >= 0.05 Sensitivity
## 6 DRIMSeq(g) 0.566 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.110 D >= 0.05 FDR
## 2 DEXSeq 0.124 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.177 D >= 0.05 FDR
## 4 RATs(t) 0.314 D >= 0.05 FDR
## 5 SUPPA2 0.654 D >= 0.05 FDR
## 6 RATs(g) 0.673 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.726 D >= 0.05 MCC
## 2 DEXSeq 0.723 D >= 0.05 MCC
## 3 DRIMSeq(t) 0.696 D >= 0.05 MCC
## 4 DRIMSeq(g) 0.672 D >= 0.05 MCC
## 5 SUPPA2 0.538 D >= 0.05 MCC
## 6 RATs(g) 0.513 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}Prefilter > 10
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Human') %>%
filter(prefilter == 10)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.608 D >= 0.2 Sensitivity
## 2 DRIMSeq(g) 0.596 D >= 0.2 Sensitivity
## 3 DRIMSeq(t) 0.589 D >= 0.2 Sensitivity
## 4 DEXSeq 0.567 D >= 0.2 Sensitivity
## 5 RATs(g) 0.563 D >= 0.2 Sensitivity
## 6 RATs(t) 0.549 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.018 D >= 0.2 FDR
## 2 RATs(g) 0.023 D >= 0.2 FDR
## 3 DEXSeq 0.032 D >= 0.2 FDR
## 4 SUPPA2 0.050 D >= 0.2 FDR
## 5 DRIMSeq(t) 0.058 D >= 0.2 FDR
## 6 DRIMSeq(g) 0.078 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.757 D >= 0.2 MCC
## 2 DEXSeq 0.738 D >= 0.2 MCC
## 3 RATs(g) 0.737 D >= 0.2 MCC
## 4 RATs(t) 0.731 D >= 0.2 MCC
## 5 DRIMSeq(g) 0.731 D >= 0.2 MCC
## 6 DRIMSeq(t) 0.726 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.787 D >= 0.1 Sensitivity
## 2 RATs(g) 0.716 D >= 0.1 Sensitivity
## 3 RATs(t) 0.692 D >= 0.1 Sensitivity
## 4 DRIMSeq(g) 0.651 D >= 0.1 Sensitivity
## 5 DRIMSeq(t) 0.641 D >= 0.1 Sensitivity
## 6 DEXSeq 0.599 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.063 D >= 0.1 FDR
## 2 RATs(t) 0.065 D >= 0.1 FDR
## 3 DRIMSeq(t) 0.115 D >= 0.1 FDR
## 4 DRIMSeq(g) 0.183 D >= 0.1 FDR
## 5 RATs(g) 0.220 D >= 0.1 FDR
## 6 SUPPA2 0.260 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.802 D >= 0.1 MCC
## 2 SUPPA2 0.759 D >= 0.1 MCC
## 3 DEXSeq 0.746 D >= 0.1 MCC
## 4 RATs(g) 0.741 D >= 0.1 MCC
## 5 DRIMSeq(t) 0.732 D >= 0.1 MCC
## 6 DRIMSeq(g) 0.717 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.878 D >= 0.05 Sensitivity
## 2 RATs(g) 0.864 D >= 0.05 Sensitivity
## 3 RATs(t) 0.778 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.656 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.649 D >= 0.05 Sensitivity
## 6 DEXSeq 0.601 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.110 D >= 0.05 FDR
## 2 DRIMSeq(t) 0.176 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.272 D >= 0.05 FDR
## 4 RATs(t) 0.310 D >= 0.05 FDR
## 5 SUPPA2 0.565 D >= 0.05 FDR
## 6 RATs(g) 0.658 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.728 D >= 0.05 MCC
## 2 DEXSeq 0.728 D >= 0.05 MCC
## 3 DRIMSeq(t) 0.707 D >= 0.05 MCC
## 4 DRIMSeq(g) 0.676 D >= 0.05 MCC
## 5 SUPPA2 0.610 D >= 0.05 MCC
## 6 RATs(g) 0.526 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}Prefilter > 50
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Human') %>%
filter(prefilter == 50)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.632 D >= 0.2 Sensitivity
## 2 DRIMSeq(t) 0.616 D >= 0.2 Sensitivity
## 3 SUPPA2 0.600 D >= 0.2 Sensitivity
## 4 RATs(g) 0.562 D >= 0.2 Sensitivity
## 5 RATs(t) 0.549 D >= 0.2 Sensitivity
## 6 DEXSeq 0.523 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.007 D >= 0.2 FDR
## 2 RATs(g) 0.009 D >= 0.2 FDR
## 3 DRIMSeq(t) 0.039 D >= 0.2 FDR
## 4 DRIMSeq(g) 0.051 D >= 0.2 FDR
## 5 DEXSeq 0.056 D >= 0.2 FDR
## 6 SUPPA2 0.059 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.763 D >= 0.2 MCC
## 2 SUPPA2 0.750 D >= 0.2 MCC
## 3 DRIMSeq(t) 0.747 D >= 0.2 MCC
## 4 RATs(g) 0.742 D >= 0.2 MCC
## 5 RATs(t) 0.735 D >= 0.2 MCC
## 6 DEXSeq 0.699 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.776 D >= 0.1 Sensitivity
## 2 RATs(g) 0.715 D >= 0.1 Sensitivity
## 3 RATs(t) 0.692 D >= 0.1 Sensitivity
## 4 DRIMSeq(g) 0.678 D >= 0.1 Sensitivity
## 5 DRIMSeq(t) 0.665 D >= 0.1 Sensitivity
## 6 DEXSeq 0.552 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.040 D >= 0.1 FDR
## 2 RATs(g) 0.081 D >= 0.1 FDR
## 3 DRIMSeq(t) 0.087 D >= 0.1 FDR
## 4 DEXSeq 0.094 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.131 D >= 0.1 FDR
## 6 SUPPA2 0.152 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.813 D >= 0.1 MCC
## 2 SUPPA2 0.810 D >= 0.1 MCC
## 3 RATs(g) 0.807 D >= 0.1 MCC
## 4 DRIMSeq(t) 0.755 D >= 0.1 MCC
## 5 DRIMSeq(g) 0.754 D >= 0.1 MCC
## 6 DEXSeq 0.704 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.902 D >= 0.05 Sensitivity
## 2 RATs(g) 0.864 D >= 0.05 Sensitivity
## 3 RATs(t) 0.777 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.683 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.671 D >= 0.05 Sensitivity
## 6 DEXSeq 0.555 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
sensh[as.character(tmp2$tool[i])] = sensh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.128 D >= 0.05 FDR
## 2 DEXSeq 0.150 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.183 D >= 0.05 FDR
## 4 RATs(t) 0.314 D >= 0.05 FDR
## 5 SUPPA2 0.566 D >= 0.05 FDR
## 6 RATs(g) 0.673 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
fdrh[as.character(tmp2$tool[i])] = fdrh[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.738 D >= 0.05 MCC
## 2 DRIMSeq(g) 0.732 D >= 0.05 MCC
## 3 RATs(t) 0.726 D >= 0.05 MCC
## 4 DEXSeq 0.683 D >= 0.05 MCC
## 5 SUPPA2 0.623 D >= 0.05 MCC
## 6 RATs(g) 0.513 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyh[as.character(tmp2$tool[i])] = tallyh[as.character(tmp2$tool[i])] + i
mcch[as.character(tmp2$tool[i])] = mcch[as.character(tmp2$tool[i])] + i
}Cumulative rank
print( data.frame(Overall = tallyh,
Sensitivity = sensh,
FDR = fdrh,
MCC = mcch) )## Overall Sensitivity FDR MCC
## RATs(t) 74 33 21 20
## RATs(g) 100 24 38 38
## DRIMSeq(t) 93 40 21 32
## DRIMSeq(g) 111 37 37 37
## SUPPA2 94 11 49 34
## DEXSeq 95 44 23 28Fruitfly
Without prefilter
tallyd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
sensd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
fdrd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
mccd = c('RATs(t)'=0, 'RATs(g)'=0, 'DRIMSeq(t)'=0, 'DRIMSeq(g)'=0, 'SUPPA2'=0, 'DEXSeq'=0)
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Fruitfly') %>%
filter(prefilter == 0)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.615 D >= 0.2 Sensitivity
## 2 DRIMSeq(g) 0.601 D >= 0.2 Sensitivity
## 3 DRIMSeq(t) 0.595 D >= 0.2 Sensitivity
## 4 DEXSeq 0.577 D >= 0.2 Sensitivity
## 5 RATs(g) 0.535 D >= 0.2 Sensitivity
## 6 RATs(t) 0.526 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.009 D >= 0.2 FDR
## 2 RATs(g) 0.015 D >= 0.2 FDR
## 3 DEXSeq 0.015 D >= 0.2 FDR
## 4 DRIMSeq(t) 0.029 D >= 0.2 FDR
## 5 DRIMSeq(g) 0.036 D >= 0.2 FDR
## 6 SUPPA2 0.128 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.743 D >= 0.2 MCC
## 2 DRIMSeq(g) 0.741 D >= 0.2 MCC
## 3 DRIMSeq(t) 0.721 D >= 0.2 MCC
## 4 RATs(g) 0.719 D >= 0.2 MCC
## 5 SUPPA2 0.718 D >= 0.2 MCC
## 6 RATs(t) 0.710 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.765 D >= 0.1 Sensitivity
## 2 RATs(g) 0.666 D >= 0.1 Sensitivity
## 3 RATs(t) 0.653 D >= 0.1 Sensitivity
## 4 DRIMSeq(t) 0.649 D >= 0.1 Sensitivity
## 5 DRIMSeq(g) 0.640 D >= 0.1 Sensitivity
## 6 DEXSeq 0.616 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.024 D >= 0.1 FDR
## 2 DEXSeq 0.033 D >= 0.1 FDR
## 3 RATs(g) 0.043 D >= 0.1 FDR
## 4 DRIMSeq(t) 0.056 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.071 D >= 0.1 FDR
## 6 SUPPA2 0.223 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.792 D >= 0.1 MCC
## 2 RATs(t) 0.788 D >= 0.1 MCC
## 3 DEXSeq 0.761 D >= 0.1 MCC
## 4 SUPPA2 0.755 D >= 0.1 MCC
## 5 DRIMSeq(g) 0.750 D >= 0.1 MCC
## 6 DRIMSeq(t) 0.743 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.801 D >= 0.05 Sensitivity
## 2 SUPPA2 0.786 D >= 0.05 Sensitivity
## 3 RATs(t) 0.754 D >= 0.05 Sensitivity
## 4 DRIMSeq(t) 0.652 D >= 0.05 Sensitivity
## 5 DRIMSeq(g) 0.642 D >= 0.05 Sensitivity
## 6 DEXSeq 0.621 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DEXSeq 0.052 D >= 0.05 FDR
## 2 DRIMSeq(t) 0.072 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.099 D >= 0.05 FDR
## 4 RATs(t) 0.118 D >= 0.05 FDR
## 5 RATs(g) 0.222 D >= 0.05 FDR
## 6 SUPPA2 0.248 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.804 D >= 0.05 MCC
## 2 RATs(g) 0.781 D >= 0.05 MCC
## 3 DEXSeq 0.756 D >= 0.05 MCC
## 4 SUPPA2 0.752 D >= 0.05 MCC
## 5 DRIMSeq(g) 0.738 D >= 0.05 MCC
## 6 DRIMSeq(t) 0.736 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}Prefilter > 10
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Fruitfly') %>%
filter(prefilter == 10)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.635 D >= 0.2 Sensitivity
## 2 SUPPA2 0.618 D >= 0.2 Sensitivity
## 3 DRIMSeq(t) 0.614 D >= 0.2 Sensitivity
## 4 DEXSeq 0.588 D >= 0.2 Sensitivity
## 5 RATs(g) 0.535 D >= 0.2 Sensitivity
## 6 RATs(t) 0.526 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.008 D >= 0.2 FDR
## 2 RATs(g) 0.009 D >= 0.2 FDR
## 3 DEXSeq 0.018 D >= 0.2 FDR
## 4 SUPPA2 0.026 D >= 0.2 FDR
## 5 DRIMSeq(t) 0.029 D >= 0.2 FDR
## 6 DRIMSeq(g) 0.033 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.762 D >= 0.2 MCC
## 2 SUPPA2 0.761 D >= 0.2 MCC
## 3 DEXSeq 0.749 D >= 0.2 MCC
## 4 DRIMSeq(t) 0.727 D >= 0.2 MCC
## 5 RATs(g) 0.722 D >= 0.2 MCC
## 6 RATs(t) 0.711 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.772 D >= 0.1 Sensitivity
## 2 DRIMSeq(g) 0.683 D >= 0.1 Sensitivity
## 3 DRIMSeq(t) 0.682 D >= 0.1 Sensitivity
## 4 RATs(g) 0.666 D >= 0.1 Sensitivity
## 5 DEXSeq 0.663 D >= 0.1 Sensitivity
## 6 RATs(t) 0.653 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.021 D >= 0.1 FDR
## 2 RATs(g) 0.036 D >= 0.1 FDR
## 3 DEXSeq 0.058 D >= 0.1 FDR
## 4 DRIMSeq(t) 0.066 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.073 D >= 0.1 FDR
## 6 SUPPA2 0.136 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.801 D >= 0.1 MCC
## 2 RATs(g) 0.795 D >= 0.1 MCC
## 3 RATs(t) 0.789 D >= 0.1 MCC
## 4 DEXSeq 0.779 D >= 0.1 MCC
## 5 DRIMSeq(g) 0.773 D >= 0.1 MCC
## 6 DRIMSeq(t) 0.753 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.801 D >= 0.05 Sensitivity
## 2 SUPPA2 0.795 D >= 0.05 Sensitivity
## 3 RATs(t) 0.754 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.685 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.684 D >= 0.05 Sensitivity
## 6 DEXSeq 0.670 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.086 D >= 0.05 FDR
## 2 DEXSeq 0.088 D >= 0.05 FDR
## 3 DRIMSeq(g) 0.098 D >= 0.05 FDR
## 4 RATs(t) 0.116 D >= 0.05 FDR
## 5 SUPPA2 0.171 D >= 0.05 FDR
## 6 RATs(g) 0.218 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.805 D >= 0.05 MCC
## 2 SUPPA2 0.794 D >= 0.05 MCC
## 3 RATs(g) 0.783 D >= 0.05 MCC
## 4 DEXSeq 0.769 D >= 0.05 MCC
## 5 DRIMSeq(g) 0.761 D >= 0.05 MCC
## 6 DRIMSeq(t) 0.742 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}Prefilter > 50
tmp <- mega %>%
filter(tool != 'DEXsgr') %>%
filter(!reproducibility %in% c('Qrep >= 0.9 Rrep >= 0.75', 'Qrep >= 0.85 Rrep >= 0.65', 'Qrep >= 0.8 Rrep >= 0.55')) %>%
filter(quantification == 'Salmon') %>%
filter(species == 'Fruitfly') %>%
filter(prefilter == 50)
# D>=0.2
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(g) 0.647 D >= 0.2 Sensitivity
## 2 DRIMSeq(t) 0.631 D >= 0.2 Sensitivity
## 3 SUPPA2 0.628 D >= 0.2 Sensitivity
## 4 DEXSeq 0.556 D >= 0.2 Sensitivity
## 5 RATs(g) 0.535 D >= 0.2 Sensitivity
## 6 RATs(t) 0.526 D >= 0.2 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(g) 0.007 D >= 0.2 FDR
## 2 RATs(t) 0.008 D >= 0.2 FDR
## 3 SUPPA2 0.012 D >= 0.2 FDR
## 4 DRIMSeq(t) 0.024 D >= 0.2 FDR
## 5 DRIMSeq(g) 0.026 D >= 0.2 FDR
## 6 DEXSeq 0.030 D >= 0.2 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.2') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.777 D >= 0.2 MCC
## 2 DRIMSeq(g) 0.766 D >= 0.2 MCC
## 3 DRIMSeq(t) 0.728 D >= 0.2 MCC
## 4 DEXSeq 0.723 D >= 0.2 MCC
## 5 RATs(g) 0.722 D >= 0.2 MCC
## 6 RATs(t) 0.711 D >= 0.2 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.1
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.781 D >= 0.1 Sensitivity
## 2 DRIMSeq(g) 0.703 D >= 0.1 Sensitivity
## 3 DRIMSeq(t) 0.700 D >= 0.1 Sensitivity
## 4 RATs(g) 0.666 D >= 0.1 Sensitivity
## 5 RATs(t) 0.653 D >= 0.1 Sensitivity
## 6 DEXSeq 0.612 D >= 0.1 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 RATs(t) 0.017 D >= 0.1 FDR
## 2 RATs(g) 0.023 D >= 0.1 FDR
## 3 DRIMSeq(t) 0.054 D >= 0.1 FDR
## 4 DEXSeq 0.058 D >= 0.1 FDR
## 5 DRIMSeq(g) 0.061 D >= 0.1 FDR
## 6 SUPPA2 0.070 D >= 0.1 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.1') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.843 D >= 0.1 MCC
## 2 RATs(g) 0.801 D >= 0.1 MCC
## 3 RATs(t) 0.791 D >= 0.1 MCC
## 4 DRIMSeq(g) 0.784 D >= 0.1 MCC
## 5 DRIMSeq(t) 0.757 D >= 0.1 MCC
## 6 DEXSeq 0.747 D >= 0.1 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}
# D>=0.05
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'Sensitivity') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.817 D >= 0.05 Sensitivity
## 2 RATs(g) 0.801 D >= 0.05 Sensitivity
## 3 RATs(t) 0.754 D >= 0.05 Sensitivity
## 4 DRIMSeq(g) 0.704 D >= 0.05 Sensitivity
## 5 DRIMSeq(t) 0.702 D >= 0.05 Sensitivity
## 6 DEXSeq 0.616 D >= 0.05 Sensitivityfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
sensd[as.character(tmp2$tool[i])] = sensd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'FDR') %>%
arrange(value) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 DRIMSeq(t) 0.063 D >= 0.05 FDR
## 2 DRIMSeq(g) 0.076 D >= 0.05 FDR
## 3 DEXSeq 0.083 D >= 0.05 FDR
## 4 RATs(t) 0.118 D >= 0.05 FDR
## 5 SUPPA2 0.149 D >= 0.05 FDR
## 6 RATs(g) 0.222 D >= 0.05 FDRfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
fdrd[as.character(tmp2$tool[i])] = fdrd[as.character(tmp2$tool[i])] + i
}
tmp2 <- tmp %>%
filter(effect == 'D >= 0.05') %>%
filter(type == 'MCC') %>%
arrange(desc(value)) %>%
dplyr::select(tool, value, effect, type)
print(tmp2)## tool value effect type
## 1 SUPPA2 0.822 D >= 0.05 MCC
## 2 RATs(t) 0.804 D >= 0.05 MCC
## 3 RATs(g) 0.781 D >= 0.05 MCC
## 4 DRIMSeq(g) 0.777 D >= 0.05 MCC
## 5 DRIMSeq(t) 0.753 D >= 0.05 MCC
## 6 DEXSeq 0.739 D >= 0.05 MCCfor (i in 1:nrow(tmp2)) {
tallyd[as.character(tmp2$tool[i])] = tallyd[as.character(tmp2$tool[i])] + i
mccd[as.character(tmp2$tool[i])] = mccd[as.character(tmp2$tool[i])] + i
}Cumulative rank
Cumulative rank
print( data.frame(Overall = tallyd,
Sensitivity = sensd,
FDR = fdrd,
MCC = mccd) )## Overall Sensitivity FDR MCC
## RATs(t) 90 41 19 30
## RATs(g) 85 29 29 27
## DRIMSeq(t) 104 32 28 44
## DRIMSeq(g) 98 26 39 33
## SUPPA2 82 14 47 21
## DEXSeq 108 47 27 34Both datasets
tmp <- data.frame(Cumulative_Rank= c(tallyh + tallyd, sensh + sensd, fdrh + fdrd, mcch + mccd),
Type=factor(rep(c('Overall', 'Sensitivity', 'FDR', 'MCC'), each=6), levels=c('Overall', 'Sensitivity', 'FDR', 'MCC')),
Tool=c(names(tallyh), names(sensh), names(fdrh), names(mcch)) )
p <-ggplot(tmp, aes(x=Tool, fill=Type, y=Cumulative_Rank)) + gth +
geom_bar(stat='identity', position='dodge') +
scale_y_reverse() +
scale_x_discrete(position = "top") +
scale_fill_manual(values=c(Overall='grey70', Sensitivity='orange', FDR='red', MCC='blue'))
print(p)